c/* @(#)optp.f 1.1 latest revision 6/23/89 14:55:12 */ subroutine OPTP ( kpts, n1, lpia, lmaxin, lborn, nm1, nm2, aovera) ************************************************************************ c *** May 82 version with ReU and ImU separate... c *** see dimensions + equivalence statements c *** Also, see R.H. Landau's program LPOTT, 1981 c Changed by RHL and LB 28 Nov 88 for 32 grid points and 100 partial waves c Changed by LB 10 Dec 88 - ur(33,33,6) -> ur(49,49,6) c ui(33,33,6) -> ui(49,49,6) c after run with 48 pts and 48 pw c *** Need to reset dimensions and equivalences as max # grid increase c should use n = ngp + 1 for scattering. c ff(ndimf/6, 6) also needed ,and equivalenced to f(1) c For equivalence need length + 1 to overlay c Also see initialize statements which now change with dimensions c *** For spin 0 * 0.5 calculate the partial wave decomposition c of rho (the form factor). c rho(nff,i) is the i-1 th partial wave amplitude of the form factor c for a value of k and kp. c *** nff determines which f.f. is used c nff = 1 proton matter = 2 neutron matter c = 3 proton spin = 4 neutron spin c *** REVISIONS *** c c MS (12/4/87) c Changed nuc0 to a logical type variable called spin0 c c *** Do all l values on 1st call, then store on tape9, ft is for temp c implicit real*8 (a-h, k, o-z) real*8 imtij, Mp, MN, k, kp, kpts, k1 Logical spin0 c c %%% changed dimension of f(600,6) to f(3201,6) to match fr(98,98),fi(98,98) c %%% in common /fopt/ with dummy(2) added - LB 10 Dec 88 dimension Vcoull(100) dimension nifty(20) dimension bnuc(40), bnucf(40), bn(40), bnf(40) dimension retij(40), imtij(40) dimension kpts(*) dimension Ur(97,97,8), Ui(97,97,8) dimension Vr(7), Vi(7), f(100,100,220,7), Vll(100,2,7) dimension dummy(2) common /bandt/ bnuc, bn, bnucf, bnf, retij, imtij common /params/ hbarc, pi, Mp, MN, nz, na, nes, nwaves common /switch/ nifty common /nlspfl/ Ur, Ui common /parallel/ f c common /foptp/ f, dummy common /sizes/ achp, acmp, wsp, achn, acmn, wsn common /ranges/ Rcoul, Rcut c equivalence (f(1,1), Vr(1)) c equivalence (f(7,1), Vi(1)), (f(1,1), ff(1,1)) data Vcoull /100 * 0.0d+00/ data nota /1/ data nint /0/ data maxLs / 100 / c >>> FIRST EXECUTABLE STATEMENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< IF (lpia .ne. 1) go to 250 c *** 1st call to optp, calculate all U's c *** Initialize (can't use data for this due to commons) c %%% changed dimension of f(600,6) to f(3201,6) to match fr(98,98),fi(98,98) c %%% so that initialization of f must be changed nint = 1 i1 = nm1 * nm1 * 4 c IF (i1 .gt. 3201) i1 = 3201 c do 7 i = 1,6 c do 7 j = 1,i1 c f(j,i) = 0. c 7 continue k = -200. spin0 = .False. c *** Special 0 * 0.5 spin case, use old t * rholl (the partial wave c *** projection of the form factor). IF (nifty(6) .eq. 3 .OR. nifty(6) .eq. 0) then spin0 = .True. ENDIF call TNCM( k ) IF (kpts(n1) .gt. 0.) then k1 = kpts(n1) else k1 = kpts(1) ENDIF write(6,410) k1 n2 = 2 * n1 nn = na - nz c *** Make sure lmax = 1 case does not run into indef on initial runs ldmax = lmaxin IF (ldmax .eq. 1) ldmax = 2 rhl3 = (achp/hbarc)**2/2. rhl4 = (achn/hbarc)**2/2. c *** Added terms for C-13 spin...may want to change size here rhl5 = (wsn/hbarc)**2/2. as = 1./6. IF (nz .ne. 0) ap = (nz - 2.)/6./nz * (acmp/achp)**2 IF (nz .le. 2) ap = 0. IF (nn .ne. 0) an = (nn - 2.)/6./nn * (acmn/achn)**2 IF (nn .le. 2) an = 0. c *** Special C-13, add in p1/2 n and rescale IF (na .eq. 13) an = (nn - 3.)/6./(nn - 1.) * (acmn/achn)**2 c *** Do loop over all grid points c rewind 9 IF (nifty(6) .eq. 6) nifty(17) = 1 aovera = (na - 1.)/na IF (nifty(17) .eq. 1) aovera = 1. IF (na .le. 2) aovera = 1. c *** Set U(strong) = 0 IF (nifty(6) .eq. 4) aovera = 0. xhl = aovera * 2. * pi * pi aparam = 1. IF (nifty(1) .eq. 9) aparam = 0. IF (Rcoul .lt. 1.e-10 ) Rcoul = 1.e-6 aparam = (1.5 * nz * aparam/Rcoul) * hbarc * 7.29735e-3 xhl1 = Rcoul/hbarc xhl2 = Rcut/hbarc do 240 i1 = 1,n1 do 240 i2 = 1,n1 k = kpts(i2) kp = kpts(i1) lmax = ldmax + 4 c *** New version....outer do loop over kp, k; inner over all l's xMN = MN c *** Do loop over l *** do 230 ldum = 1,ldmax IF ( (ldum .gt. lborn) $ .AND. $ ( (i1 + i2) .ne. (2 * n1) )) go to 230 j1 = ldum Urc = 0. Uic = 0. Urs = 0. Uis = 0. c *** Set up potential needed for this l in integral equation c *** N.B. the f is + ImV as stored. npot1 = 2 * ldum IF (ldmax .gt. maxLs) write(6,998) c *** c *** Vll(neles, re/im, nspin-mp) = Vll(50,2,6) c *** lmas2 = ldmax + 2 IF (ldum .eq. 1) then if ( spin0 ) then call Vopt0A( Vll, k, kp, lmas2, k1 ) else call VStoVL( Vll, k, kp, lmas2, k1 ) endif endif c *** Proton c *** c *** Nifty(6) = 8, full 1/2 * 1/2 c *** = 3, 0 * 1/2 (can't yet do with Paez) c *** see codes latter c *** c *** This is for singlet ld2 = ldum + 2 if ( spin0 ) then f(i1,i2,npot1-1,1) = Vll( ldum, 1, 1 ) f(i1,i2,npot1,1) = Vll( ldum, 2, 1 ) else f(i1,i2,npot1-1,1) = Vll( ldum, 1, 6 ) f(i1,i2,npot1,1) = Vll( ldum, 2, 6 ) endif ff1 = f(i1,i2,npot1-1,1) ff2 = f(i1,i2,npot1,1) IF (nifty(6) .ne. 3 $ .AND. $ nifty(6) .ne. 8) go to 226 c *** c *** Coupled channels pp c *** Store matrices with new, nonlinear forms c *** c *** RHL code for proton He3 nspin c *** 1- singlet 2- triplet l=j 3- triplet j=l+1 (l,l) c *** 4- j=l+1 (l+2,l) 5- j=l+1 (l+2,l+2) 6- j=l+1 (l,l+2) c *** N.B. switch l values so nspin = 3, 4, 5, 6 all same j - as in main c *** M.P. code versus RHL 1=6 2=4 3=3 4=2 5=5 6=1 c *** c********* tvm ST Mixing 9/8/92 Yet another convention. ************* c 7-triplet l = j 2- triplet singlet 1s 8-singlet triplet s1 *********************************************************************** c 2's to sevens in f( , )'s below c 4's and 6's interchanges in code below in f( , )'s if ( spin0 ) then f(i1,i2,npot1-1,2) = Vll( ldum, 1, 2 ) f(i1,i2,npot1,2) = Vll( ldum, 2, 2 ) else c twelve ff's to f's tm 9/24/90 f(i1,i2,npot1-1,7) = Vll(ldum,1,4) f(i1,i2,npot1,7) = Vll(ldum,2,4) c *** Special Paez coding store r0(11) in r0(00) IF (ldum.eq.1) then f(i1,i2,npot1-1,7) = Vll(2,1,5) f(i1,i2,npot1,7) = Vll(2,2,5) endif f(i1,i2,npot1-1,4) = Vll(ld2,1,1) f(i1,i2,npot1,4) = Vll(ld2,2,1) f(i1,i2,npot1-1,3) = Vll(ldum,1,3) f(i1,i2,npot1,3) = Vll(ldum,2,3) f(i1,i2,npot1-1,5) = Vll(ld2,1,5) f(i1,i2,npot1,5) = Vll(ld2,2,5) f(i1,i2,npot1-1,6) = Vll(ldum,1,2) f(i1,i2,npot1,6) = Vll(ldum,2,2) c *********** New TS mixing tvm 9/8/92 *************************** f(i1,i2,npot1-1,2) = Vll(ldum,1,7) f(i1,i2,npot1,2) = Vll(ldum,2,7) endif c *** Not special Paez coding any longer 226 nnn = npot1 - 1 IF ( kp .eq. k1 $ .AND. $ kp .eq. k $ .AND. $ ldum .eq. 1) write(6,410) k1 IF (kp .ne. k1 .OR. ldum .ne. 1) go to 230 c *** Write out potential matrices on printer (half off-shell) c ******** tm 10/05/90 rm write statement, output too big c write(6,400) k, ff1, ff2 n = nifty(6) IF (n .eq. 3 .OR. n .eq. 8) then c tm 9/24/90 ff's to f's in 2 writes c write(6,401) (i, f(nnn,i), f(npot1,i), i=2,2,1) endif IF (n .eq. 8) then c write(6,401) (i, f(nnn,i), f(npot1,i), i=3,6,1) endif c *** Do loop over ldum ENDS 230 continue npot1m = 2 * lborn IF (ldmax .lt. lborn) npot1m = 2 * ldmax IF ((i1 .eq. i2) .AND. (i1 .eq. n1)) npot1m = 2 * ldmax c *** New unformatted read for time saving c IF (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) then c write(9) (f(npot1,1), npot1=1,npot1m) c endif c IF (nifty(6) .eq. 3) then c write(9) ( ( f(npot2-1,i), f(npot2,i), i=1,2 ), c $ npot2 = 2,npot1m,2 ) c endif c IF (nifty(6) .eq. 8) then c tm 9/24/90 ff's to f's c write(9) ( ( f(npot2-1,i), f(npot2,i), i= 1,7 ), c $ npot2 = 2,npot1m,2 ) c endif ******** Write variables for checking ************************ do 7170 npot2 = 2, npot1m, 2 do 7171 i = 1, 7 c write(7,7172) i1,i2,f(i1,i2,npot2-1,i),f(i1,i2,npot2,i) 7171 continue 7170 continue 7172 format(2i5,2e12.4) ******************************************************************* c *** Do loop over kp and k ENDS c *** tm 9/18/90 insert write variables c *** tm 9/29/90 print diagonal c if (lpia .ne. 1) go to 235 if (i2 .ne. i1) go to 235 c write(6,1010) lid,lpia c tm 9/24/90 ff's to f's in next line c write(7,1002) kpts(i1), f(1,1), f(2,1) 235 continue if (i2.eq.n1-1) then write(8,238) i1, f(i1,5,6,6) endif 238 format(i8,e12.4) 240 continue c *** c *** Initialization over, come here on all other calls for read c *** 250 rewind 9 c *** Potentials all stored, read in the values for current lpia c *** npts = # of i1, i2 combos (followed by U's for all l's) npts = (n1 * (1 + n1))/2 IF (lpia .ne. 1) go to 260 c *** First call for read, setup parameters c *** N.B....the values of Vl(i,j) for all l but fixed c *** i, and j are on one record via unformatted write c *** need to skip to new record to get new ij, nspin = 1 npot = 2 IF (nifty(6) .eq. 8) nspin = 7 IF (nifty(6) .ne. 3) go to 260 nspin = 2 c *** Spin 0 x 1/2 case, 4 V's now adjacent for each l vs 2 for 0 x0 npot = 2 * npot 260 ldum = 2 * (lpia - 1) * nspin IF (nifty(6) .eq. 8) go to 345 IF (lpia .gt. lborn) go to 290 c *** Do loop over npts (i.e. records) c *** N.B. ldum contains both # of l's and nspin (potential types) to skip do 280 npt = 1,npts c IF (lpia .ne. 1) go to 270 c read(9,end=10000) (f(npt,i), i=1,npot) 10000 go to 280 c *** Space to l value needed c270 read(9,end=280) (rhl,i=1,ldum), (f(npt,i), i=1,npot) 280 continue go to 310 c *** lpia .gt. lborn, read only last on-shell point 290 nskip = npts - 1 do 300 i = 1,nskip read(9,end=300) 300 continue c read(9,end=310) (rhl,i=1,ldum), (f(npts,i), i=1,npot) c *** Set up potential matrix 310 istart = 1 IF (lpia .gt. lborn) istart = n1 do 340 i1 = istart,n1 do 340 i2 = istart,n1 npt = ((i1 - 1) * i1)/2 + i2 c *** RHL change FEB 82 for bs and MAY 82 for 1/2 x 1/2 c IF (nifty(6) .eq. 7) f(npt,2) = 0. c Ur(i1,i2,1) = f(npt,1) c Ui(i1,i2,1) = f(npt,2) c *** Spin-dependent potentials IF (nifty(6) .ne. 3) go to 330 c Ur(i1,i2,2) = f(npt,3) c Ui(i1,i2,2) = f(npt,4) IF (nifty(6) .eq. 7) Ui(i1,i2,2) = 0. c *** Set up potential with U(kp,k) = U(k,kp)..symmetric U 330 IF (i1 .eq. i2) go to 340 c Ur(i2,i1,1) = Ur(i1,i2,1) c Ui(i2,i1,1) = Ui(i1,i2,1) IF (nifty(6) .ne. 3) go to 340 Ui(i2,i1,2) = Ui(i1,i2,2) Ur(i2,i1,2) = Ur(i1,i2,2) 340 continue return c *** c *** Special spin 1/2 x 1/2 coupled channels c *** nspin = 1 singlet 2-triplet j=l 3/4 j=l+1 diag+couple c *** 5/6 j=l-1 diag+couple c *** npt is location in triangular matrix stored linearly c *** rearrange with (i,j) index, supermatrix formed in main 345 istart = 1 IF (lpia .le. lborn) go to 421 c *** lpia gt lborn, read only last on-shell point nskip = npts - 1 do 10017 i = 1,nskip c read(9,end=10017) 10017 continue istart = n1 c *** Regular or Born 421 do 346 i1 = istart,n1 do 346 i2 = istart,n1 npot4 = 2*lpia npot5 = npot4-1 do 346 i = 1, 7 Ur(i1,i2,i) = f(i1,i2,npot5,i) Ui(i1,i2,i) = f(i1,i2,npot4,i) IF (lpia .eq. 1) then c read(9) (Ur(i1,i2,i), Ui(i1,i2,i), i=1,7) else c read(9) (rhl,i=1,ldum), (Ur(i1,i2,i), Ui(i1,i2,i), i=1,7) endif c ********* tvm 9/8/92 TS-mixing ************************************* Ur(i1,i2,8) = Ur(i1,i2,2) Ui(i1,i2,8) = Ui(i1,i2,2) ********************************************************************** IF (i1 .eq. i2) go to 346 c *** Build in non-diagonal elements via symmetry c do 348 i = 1,8 c ili = i c IF (i .eq. 4) ili = 6 c IF (i .eq. 6) ili = 4 c Ur(i2,i1,i) = Ur(i1,i2,ili) c 348 Ui(i2,i1,i) = Ui(i1,i2,ili) 346 continue ******* Second Check point ********************************** c write(7,*) ' ' do 8181 i1 = 1, n1 do 8080 i2 = 1, n1 do 8282 i = 1, 7 c write(7,8383) i1, i2,Ur(i1,i2,i),Ui(i1,i2,i) 8282 continue 8080 continue 8181 continue 8383 format(2i5,2e12.4) c *** tm 9/19/90 print variables if (lpia .ne. 1) go to 349 c write(6,1200) lpia,lborn c write(6,1210) Ur(lid,lid,1), Ui(lid,lid,1) 349 continue return c *** c350 format(' C13 neutron spin ff for k = kp = ',e15.5/10(10x,6e15.6/)) c400 format (3x, 2e19.4, 22x, e19.4) c401 format(1x, i2, 19x, e19.4, 41x, e19.4) 410 format('0','OPTICAL POTENIAL (OPTP) for k = k0 =', $ f13.3,' MeV and l = 0', $ /, 16x, 'kp', 17x, 'ReV', 38x, '-ImV', /) 998 format(' ','ldmax too big in OPTP at 998**********') 1000 format(' ???1 ldum= ',i2,' k=',e11.4,' kp=',e11.4) c1001 format(' Vll(ldum,1,1)=',e12.6,' Vll(ldum,2,1)='e12.6) 1002 format(' ',e10.4,' ',e10.4,' ',e10.4) 1005 format(' Vll(ldum,1,6)=',e12.6,' Vll(1,2,6)=',e12.6) c1010 format(' ???2 lid=',i3,' lpia=',i3) c1015 format(' f(1,1)= ',e10.4,' f(2,1)= ',e10.4) c1200 format(' ???3 lpia= ',i3,'lborn= ',i3) c1210 format(' Ur(lid,lid,1)=',e13.7,' Ui(lid,lid,1)=',e13.7) end