C ************ program cdxs C ************ C Stefan Typel C 2001/09/21, version 3.21 C implicit none C integer*4 file_in,file_out,file_gs,file_sc,file_dbde & ,file_ph,file_cp,file_phl,file_cpl,file_sf & ,file_cd,file_cdl,file_dbdel,file_alm,file_ev & ,nf_max,l_max,nth_max,n_max,nb_max,nx_max,nres_max parameter (file_in=10,file_out=11,file_alm=12,file_ev=24 & ,file_gs=13,file_sc=14,file_dbde=15,file_dbdel=16 & ,file_ph=17,file_phl=18,file_cp=19 & ,file_cpl=20,file_sf=21,file_cd=22,file_cdl=23 & ,nf_max=20,l_max=3,n_max=5000,nb_max=5 & ,nth_max=90,nx_max=501,nres_max=5) integer*4 i1,i2,i3,itmp1,itmp2,lamax & ,nth,ne,np,nr,nc,nx,ndum,n_thbc,n_phbc,nev & ,case_m,case_f,case_r,case_cs,case_th,case_v & ,z_x,a_x,z_a,a_a,z_b,a_b,z_c,a_c,nres(nf_max) & ,n_i,nf,incl(nf_max),el(0:l_max),lambda_f(nf_max) & ,l_i,l_f(nf_max),lmin,lmax,tja,tjmin,tjmax,ts real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu,h2m & ,m_x,m_a,m_b,m_c,mu_ax,mu_bc,zze2,xi_min & ,j_a,j_b,j_c,s,j_f(nf_max),sc(l_max),phase(l_max) & ,e_i,e_a,er_min,er_max,p_min,p_max,t_p,de,e_x & ,h,r_max,r0,r_x,r_a,a,v,gamma,theta & ,th_min,th_max,dth,b_gr,th_gr,b_min & ,th_bc_min,th_bc_max,ph_bc_min,ph_bc_max & ,thbcv(0:nth_max),phbcv(0:nth_max) & ,dth_bc,th_bc,dph_bc,ph_bc,p_bc,dofs,dsdo,d3s & ,vc_i,rc_i,ac_i,fso_i,rso_i,aso_i,rrc_i,rrso_i & ,vc_f(nf_max),rc_f(nf_max),ac_f(nf_max) & ,fso_f(nf_max),rso_f(nf_max),aso_f(nf_max) & ,rrc_f(nf_max),rrso_f(nf_max),rv(l_max,0:n_max) & ,thv(0:nth_max),wf_i(nb_max,0:n_max),z_eff(l_max) & ,eb(nb_max),delta(nf_max,0:nx_max) & ,xv(0:nx_max),wf_f(0:n_max),ddelta(nf_max,0:nx_max) & ,dp,xp(-nx_max:nx_max),dsdp(-nx_max:nx_max) & ,tmp1,tmp2,tmp3,tmp4,hi1,hi2,h2,twh,twh2,prec,acc & ,dbde(0:nx_max,nf_max),dbdel(0:nx_max,l_max) & ,sig_ph(0:nx_max,nf_max) & ,sig_cp(0:nx_max,nf_max) & ,sig_phl(0:nx_max,0:l_max),sig_cpl(0:nx_max,0:l_max) & ,sig_pht(0:nx_max),sig_cpt(0:nx_max) & ,sfac(0:nx_max),sfacl(0:nx_max,0:l_max) & ,dsde(0:nx_max,nf_max) & ,dsdel(0:nx_max,0:l_max),dsdet(0:nx_max) & ,dsdt(0:nth_max,nf_max) & ,dsdtl(0:nth_max,0:l_max),dsdtt(0:nth_max) & ,sig(nf_max),sigl(0:l_max),sigt & ,d2sde(0:nx_max,nf_max) & ,d2sdel(0:nx_max,0:l_max),d2sdet(0:nx_max) & ,d2sdt(0:nth_max,nf_max) & ,d2sdtl(0:nth_max,0:l_max),d2sdtt(0:nth_max) & ,sigtp(0:nth_max,0:nth_max) & ,eres(nf_max,nres_max),gres(nf_max,nres_max) complex*16 i,wffac,idum,idum2,z_effs(l_max) character*3 cel(3) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /wave_i/ eb,wf_i common /wave_f/ xv,wf_f,wffac common /fstates/ z_effs,j_f,l_f,lambda_f,incl,el common /pot_f/ vc_f,rrc_f,ac_f,fso_f,rrso_f,aso_f common /par_bc/ mu_bc,h2m,zze2 common /par_ax/ z_x,z_a,m_x,m_a,mu_ax,t_p common /phases/ delta,ddelta common /redtrans/ dbde,dbdel common /sigma2/ sig_ph,sig_cp,sig_phl,sig_cpl & ,sig_pht,sig_cpt,sfac,sfacl common /angles/ thv common /angbc/ thbcv,phbcv common /sig_cd/ sig,sigl,sigt common /ds_cde/ dsde,dsdel,dsdet common /ds_cdt/ dsdt,dsdtl,dsdtt common /d2s_cde/ d2sde,d2sdel,d2sdet common /d2s_cdt/ d2sdt,d2sdtl,d2sdtt common /d3s_cd/ sigtp common /reson/ eres,gres,nres common /momenta/ xp common /momdistr/ dsdp C cel(1) = 'E 1' cel(2) = 'E 2' cel(3) = 'E 3' C C ***** begin read input ***** C open(file_in,status='unknown') open(file_out,status='unknown') C read(file_in,1007) read(file_in,*) z_x read(file_in,*) a_x read(file_in,*) m_x read(file_in,*) read(file_in,*) z_a read(file_in,*) a_a read(file_in,*) m_a read(file_in,*) j_a read(file_in,*) read(file_in,*) z_b read(file_in,*) a_b read(file_in,*) m_b read(file_in,*) j_b read(file_in,*) read(file_in,*) z_c read(file_in,*) a_c read(file_in,*) m_c read(file_in,*) j_c read(file_in,*) read(file_in,*) r0 read(file_in,1001) read(file_in,*) case_m read(file_in,1001) read(file_in,*) s read(file_in,1001) read(file_in,*) h read(file_in,*) r_max read(file_in,1002) read(file_in,*) n_i read(file_in,*) read(file_in,*) l_i read(file_in,*) e_i read(file_in,*) case_v read(file_in,1010) read(file_in,*) vc_i,rc_i,ac_i,fso_i,rso_i,aso_i read(file_in,1002) read(file_in,*) case_f read(file_in,1001) read(file_in,*) nf if (nf.gt.nf_max) call err_msg(file_out,1) read(file_in,1016) do 100 i1=1,nf,1 read(file_in,*) incl(i1),lambda_f(i1),l_f(i1),j_f(i1) & ,vc_f(i1),rc_f(i1),ac_f(i1) & ,fso_f(i1),rso_f(i1),aso_f(i1) if (lambda_f(i1).gt.l_max) lambda_f(i1)=0 100 continue read(file_in,1002) read(file_in,*) el(0) = 0 read(file_in,*) (el(i1),i1=1,l_max) lamax = 0 do 105 i1=1,l_max,1 if (el(i1).eq.1) lamax=i1 105 continue read(file_in,*) (sc(i1),i1=1,l_max) read(file_in,*) (phase(i1),i1=1,l_max) read(file_in,1002) read(file_in,*) e_a read(file_in,1001) read(file_in,*) case_r read(file_in,1003) read(file_in,*) case_cs read(file_in,1028) read(file_in,*) case_th read(file_in,*) th_min read(file_in,*) th_max read(file_in,*) nth if (nth.gt.nth_max) call err_msg(file_out,2) read(file_in,*) read(file_in,*) er_min read(file_in,*) er_max read(file_in,*) ne if (ne.gt.nx_max) call err_msg(file_out,3) read(file_in,*) read(file_in,*) th_bc_min read(file_in,*) th_bc_max read(file_in,*) n_thbc read(file_in,*) read(file_in,*) ph_bc_min read(file_in,*) ph_bc_max read(file_in,*) n_phbc read(file_in,*) read(file_in,*) p_max read(file_in,*) np if ((2*np+1).gt.nx_max) call err_msg(file_out,4) read(file_in,*) read(file_in,*) nev close(file_in) C C ***** end read input ***** C C +++++++++++++++ call init_const C +++++++++++++++ C if (case_cs.eq.20) then open(file_alm,status='unknown') open(file_ev,status='unknown') C ++++++++++++++++++++++++++++++++++++ call simul(nev,file_alm,file_ev,acc) C ++++++++++++++++++++++++++++++++++++ close(file_alm) close(file_ev) goto 999 endif C c nth = nth+mod(nth,2) c ne = ne+mod(ne,2) c np = np+mod(np,2) ndum = nth_max if (ndum.gt.n_max) call err_msg(file_out,5) C nr = idnint(r_max/h) c nr = nr+mod(nr,2) if (nr.gt.n_max) call err_msg(file_out,8) r_max = h*dfloat(nr) C C ***** default values ***** C prec = 1.d-08 C if (m_x.le.0.d00) m_x = dfloat(a_x) if (m_x.le.0.d00) m_a = dfloat(a_a) if (m_x.le.0.d00) m_b = dfloat(a_b) if (m_x.le.0.d00) m_c = dfloat(a_c) C if (rc_i.le.0.d00) rc_i = 1.2d00 if (ac_i.le.0.d00) ac_i = 0.65d00 if (fso_i.lt.0.d00) fso_i = 0.358d00 if (rso_i.le.0.d00) rso_i = 1.2d00 if (aso_i.le.0.d00) aso_i = 0.65d00 C if (case_m.eq.1) then s = j_b j_c = 0.d00 endif C if (case_f.eq.0) then lmin = l_i lmax = l_i tja = idnint(2.d00*j_a) tjmin = tja tjmax = tja ts = idnint(2.d00*s) do 130 i1=1,l_max,1 if (el(i1).eq.1) then if (lmin.gt.iabs(l_i-i1)) lmin=iabs(l_i-i1) if (lmax.lt.iabs(l_i+i1)) lmax=iabs(l_i+i1) if (tjmin.gt.iabs(tja-2*i1)) tjmin=iabs(tja-2*i1) if (tjmax.lt.iabs(tja+2*i1)) tjmax=iabs(tja+2*i1) endif 130 continue nf = 0 do 140 i1=lmin,lmax,1 do 150 i2=tjmin,tjmax,2 if ((iabs(i2-2*i1).le.ts).and.(iabs(i2+2*i1).ge.ts)) then do 155 i3=1,l_max,1 if (el(i3).eq.1) then if ((iabs(i1-l_i).le.i3).and.(iabs(i1+l_i).le.i3).and. & (iabs(i2-2*i3).le.tja).and.(iabs(i2+2*i3).le.tja)) then nf = nf+1 incl(nf) = 1 lambda_f(nf) = i3 l_f(nf) = i1 j_f(nf) = dfloat(i2)/2.d00 vc_f(nf) = 0.d00 rc_f(nf) = 1.2d00 ac_f(nf) = 0.65d00 fso_f(nf) = 0.d00 rso_f(nf) = 1.2d00 aso_f(nf) = 0.65d00 endif endif 155 continue endif 150 continue 140 continue else do 120 i1=1,nf,1 if (incl(i1).eq.1) then if (vc_f(i1).le.0.d00) vc_f(i1) = 0.d00 if (rc_f(i1).le.0.d00) rc_f(i1) = 1.2d00 if (ac_f(i1).le.0.d00) ac_f(i1) = 0.65d00 if (fso_f(i1).lt.0.d00) fso_f(i1) = 0.d00 if (rso_f(i1).le.0.d00) rso_f(i1) = 1.2d00 if (aso_f(i1).le.0.d00) aso_f(i1) = 0.65d00 endif 120 continue endif C C ***** masses in MeV/c^2 ***** m_x = m_x*amu m_a = m_a*amu m_b = m_b*amu m_c = m_c*amu mu_ax = m_a*m_x/(m_a+m_x) mu_bc = m_b*m_c/(m_b+m_c) C C ***** parameters for system b+c ***** h2m = hbarc*hbarc/(2.d00*mu_bc) zze2 = dfloat(z_b*z_c)*e2 C C ***** total projectile energy in MeV ***** t_p = dfloat(a_a)*e_a C C ***** radii in fm ***** r_x = r0*dfloat(a_x)**0.333333333333d00 tmp1 = dfloat(a_a)**0.333333333333d00 r_a = r0*tmp1 rrc_i = rc_i*tmp1 rrso_i = rso_i*tmp1 fso_i = fso_i*rso_i do 160 i1=1,nf,1 rrc_f(i1) = rc_f(i1)*tmp1 rrso_f(i1) = rso_f(i1)*tmp1 fso_f(i1) = fso_f(i1)*rso_f(i1) 160 continue C C ***** init radial mesh ***** h2 = h*h hi1 = 12.d00/h2 hi2 = 2.d00*hi1 twh = 12.d00*h twh2 = twh*h do 170 i1=0,nr,1 rv(1,i1) = h*dfloat(i1) do 180 i2=2,l_max,1 rv(i2,i1) = rv(i2-1,i1)*rv(1,i1) 180 continue 170 continue C C ***** step 1: calculation of bound states ***** C +++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_wf_bound(nr,z_b,z_c,n_i,l_i,s,j_a,prec & ,rrc_i,ac_i,vc_i,rso_i,aso_i,fso_i,e_i,case_v,nc) C +++++++++++++++++++++++++++++++++++++++++++++++++++++ if (nc.gt.0) call err_msg(file_out,9) if (nc.lt.0) call err_msg(file_out,10) if (e_i.lt.0.d00) call err_msg(file_out,6) C C ***** step 2: calculation of scattering states ***** C ***** and reduced matrix elements ***** C C ***** effective charge numbers ***** tmp1 = m_c/(m_b+m_c) tmp2 = -m_b/(m_b+m_c) tmp3 = 1.d00 tmp4 = 1.d00 do 200 i1=1,l_max,1 tmp3 = tmp3*tmp1 tmp4 = tmp4*tmp2 c if (el(i1).eq.1) then z_eff(i1) = dfloat(z_b)*tmp3+dfloat(z_c)*tmp4 z_effs(i1) = z_eff(i1)*sc(i1)*cdexp(i*phase(i1)*radfac) c else c z_eff(i1) = 0.d00 c endif 200 continue C C ++++++++++++++++++++++++++++++ call get_coeff_d(nf,l_i,s,j_a) C ++++++++++++++++++++++++++++++ C C ***** cm energies E_bc***** h2m = hbarc*hbarc/(2.d00*mu_bc) if (case_cs.le.16) then if ((case_cs.eq.5).or.(case_cs.eq.7).or.(case_cs.eq.8).or. & (case_cs.ge.13)) then nx = 0 de = 0.d00 xv(0) = er_min er_max = er_min else if (ne.gt.0) then de = (er_max-er_min)/dfloat(ne) else de = 0.d00 er_max = er_min endif nx = ne do 210 i1=0,nx,1 xv(i1) = er_min+de*dfloat(i1) 210 continue endif else if (case_cs.eq.19) then if (ne.gt.0) then de = (er_max-er_min)/dfloat(ne) else de = 0.d00 er_max = er_min endif nx = ne do 215 i1=0,nx,1 xv(i1) = er_min+de*dfloat(i1) 215 continue else nx = np dp = p_max/dfloat(np) do 220 i1=(-np-1),(np+1),1 xp(i1) = dp*dfloat(i1) 220 continue do 225 i1=0,(np+1),1 xv(i1) = xp(i1)*xp(i1)/(2.d00*mu_bc) 225 continue er_max = xv(np) endif endif C C ++++++++++++++++++++++++++++++++++ call get_redmat(n_i,nf,nx,nr,s,de) C ++++++++++++++++++++++++++++++++++ C if (case_cs.le.3) then C +++++++++++++++++++++++++++++ call get_resonances(de,nf,nx) C +++++++++++++++++++++++++++++ endif C C ++++++++++++++++++++++++++++++ call get_dbde(nf,nx,mu_bc,j_a) C ++++++++++++++++++++++++++++++ C C ***** step 3: calculation of two-body cross sections ***** C C ++++++++++++++++++++++++++++++++++++++ call get_sigma2(nf,nx,e_i,j_a,j_b,j_c) C ++++++++++++++++++++++++++++++++++++++ C c goto 999 C C ***** step 4: calculation of three-body cross sections ***** C if (case_r.eq.0) then C ++++++++++++++++ call init_ylmpih C ++++++++++++++++ endif C ***** scattering angles ***** C ++++++++++++++++++++++++++++++ call kinematics(e_i,a,v,gamma) C ++++++++++++++++++++++++++++++ b_gr = r_x+r_a th_gr = 2.d00*datan(a/b_gr)/radfac xi_min = e_i/(hbarc*v*gamma)*(b_gr+pih*a/gamma) C if ((case_cs.eq.4).or.(case_cs.eq.6).or.(case_cs.eq.8).or. & (case_cs.eq.11).or.(case_cs.eq.12).or.(case_cs.eq.18)) then case_th = 1 nth = 0 dth = 0.d00 thv(0) = th_min th_max = th_min else if (case_th.eq.0) then th_min = 0.d00 th_max = th_gr endif if (nth.gt.0) then dth = (th_max-th_min)/dfloat(nth) else dth = 0.d00 th_max = th_min endif do 300 i1=0,nth,1 thv(i1) = th_min+dth*dfloat(i1) 300 continue endif b_min = a/dtan(th_max*radfac/2.d00) C C if (case_cs.eq.1) then C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dsdet(e_i,de,dth,nf,nx,nth,z_x,case_r,1,lamax) C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif C if (case_cs.eq.2) then C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dsdet(e_i,de,dth,nf,nx,nth,z_x,case_r,2,lamax) C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif C if (case_cs.eq.4) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d2sdet(e_i,de,dth,nf,nx,nth,z_x,case_r,1,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif C if (case_cs.eq.5) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d2sdet(e_i,de,dth,nf,nx,nth,z_x,case_r,2,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif C if ((case_cs.eq.3).or.((case_cs.ge.6).and.(case_cs.le.16))) then C ++++++++++++++++++++++++++++++ call get_coeff_b(nf,l_i,s,j_a) C ++++++++++++++++++++++++++++++ C if ((case_cs.eq.9).or.(case_cs.eq.11).or.(case_cs.eq.13).or. & (case_cs.eq.15)) then n_phbc = 0 endif if ((case_cs.eq.10).or.(case_cs.eq.12).or.(case_cs.eq.14).or. & (case_cs.eq.16)) then n_thbc = 0 endif C if (n_thbc.gt.0) then dth_bc = (th_bc_max-th_bc_min)/dfloat(n_thbc) else dth_bc = 0.d00 endif if (n_phbc.gt.0) then dph_bc = (ph_bc_max-ph_bc_min)/dfloat(n_phbc) else dph_bc = 0.d00 endif do 330 i1=0,n_thbc,1 thbcv(i1) = th_bc_min+dth_bc*dfloat(i1) 330 continue do 340 i2=0,n_phbc,1 phbcv(i2) = ph_bc_min+dph_bc*dfloat(i2) 340 continue C if (case_cs.eq.3) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d3scd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,0,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.6) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d3scd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,1,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.7) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d3scd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,2,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.8) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d3scd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,3,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.9) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,0,0,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.10) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,0,1,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.11) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,1,0,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.12) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,1,1,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.13) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,2,0,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.14) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,2,1,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.15) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,3,0,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.16) then C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,3,1,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif endif C if ((case_cs.eq.17).or.(case_cs.eq.18)) then C ++++++++++++++++++++++++++++++ call get_coeff_b(nf,l_i,s,j_a) C ++++++++++++++++++++++++++++++ if (case_cs.eq.17) then C +++++++++++++++++++++++++++++++++++++++++++ call get_dsdp(e_i,de,dp,dth & ,z_x,case_r,nf,nx,np,nth,1,lamax) C +++++++++++++++++++++++++++++++++++++++++++ endif if (case_cs.eq.18) then C +++++++++++++++++++++++++++++++++++++++++++ call get_dsdp(e_i,de,dp,dth & ,z_x,case_r,nf,nx,np,nth,2,lamax) C +++++++++++++++++++++++++++++++++++++++++++ endif endif C if (case_cs.eq.19) then C open(file_alm,status='unknown') C write(file_alm,*) m_x,m_a write(file_alm,*) m_b,m_c write(file_alm,*) e_i,t_p write(file_alm,*) th_min,th_max write(file_alm,*) er_min,er_max write(file_alm,*) nth,nx C C ++++++++++++++++++++++++++++++ call get_coeff_b(nf,l_i,s,j_a) C ++++++++++++++++++++++++++++++ C C +++++++++++++++++++++++++++++++++++++++++ call init_simul(e_i,de,dth,z_x,case_r & ,nf,nx,nth,file_alm,lamax) C +++++++++++++++++++++++++++++++++++++++++ C close(file_alm) C endif C 999 continue C ************* output ************* C open(file_cdl,status='unknown') open(file_cd,status='unknown') C write(file_out,3020) '********************************' write(file_out,3020) '********************************' write(file_out,3020) '***** CDXS 3.21 *****' write(file_out,3020) '********************************' write(file_out,3020) '***** Coulomb dissociation *****' write(file_out,3020) '***** cross sections *****' write(file_out,3020) '********************************' write(file_out,3020) '***** 2001/09/21 *****' write(file_out,3020) '********************************' write(file_out,3020) '********************************' C if (case_cs.ge.19) then C write(file_cd,*) write(file_cdl,*) endif C if (case_cs.eq.20) then write(file_out,*) write(file_out,*) 'number of events',nev write(file_out,*) & 'acceptance ratio in Metropolis algorithm',acc write(file_out,*) stop endif C write(file_out,*) write(file_out,3013) '***** parameters of nuclei *****' write(file_out,*) c write(file_out,3014) 'target','Z =',z_x,'A =',a_x write(file_out,3070) 'target X ' write(file_out,3071) 'Z =',z_x,'A =',a_x c write(file_out,3015) (m_x/amu),'u',m_x,'MeV' write(file_out,3072) 'm =',(m_x/amu),'u =',m_x,'MeV' write(file_out,*) c write(file_out,3016) 'projectile','Z =',z_a,'A =',a_a write(file_out,3070) 'projectile a' write(file_out,3071) 'Z =',z_a,'A =',a_a,'J =',j_a c write(file_out,3015) (m_a/amu),'u',m_a,'MeV' write(file_out,3072) 'm =',(m_a/amu),'u =',m_a,'MeV' write(file_out,*) c write(file_out,3016) 'fragment b','Z =',z_b,'A =',a_b write(file_out,3070) 'fragment b ' write(file_out,3071) 'Z =',z_b,'A =',a_b,'J =',j_b c write(file_out,3015) (m_b/amu),'u',m_b,'MeV' write(file_out,3072) 'm =',(m_b/amu),'u =',m_b,'MeV' write(file_out,*) c write(file_out,3016) 'fragment c','Z =',z_c,'A =',a_c write(file_out,3070) 'fragment c ' write(file_out,3071) 'Z =',z_c,'A =',a_c,'J =',j_c c write(file_out,3015) (m_c/amu),'u',m_c,'MeV' write(file_out,3072) 'm =',(m_c/amu),'u =',m_c,'MeV' C write(file_out,*) if (case_m.eq.1) then write(file_out,3073) 'single particle model' else write(file_out,3073) 'cluster model ' endif C write(file_out,*) write(file_out,3074) 'step size for discretization',h,'fm' write(file_out,3075) 'maximum radius',r_max,'fm' write(file_out,3076) 'number of points',nr C write(file_out,*) write(file_out,3002) '***** ground state of projectile *****' write(file_out,*) if (case_m.ne.1) then write(file_out,3003) 'n =',n_i,'l =',l_i,'j =',j_a,'s =',s else write(file_out,3003) 'n =',n_i,'l =',l_i,'j =',j_a endif write(file_out,*) write(file_out,3004) 'binding energy',(-e_i),'MeV' write(file_out,*) write(file_out,3005) 'potential parameters' write(file_out,3006) 'V_c','r_c','a_c','f_so','r_so','a_so' write(file_out,3007) '[MeV]','[fm]','[fm]','[fm]','[fm]' write(file_out,3008) vc_i,rc_i,ac_i,fso_i/rso_i,rso_i,aso_i C write(file_out,*) write(file_out,3009) '***** continuum states *****' write(file_out,*) write(file_out,3005) 'potential parameters' write(file_out,3010) 'c','lambda','l','j','V_c','r_c','a_c' & ,'f_so','r_so','a_so' write(file_out,3012) '[MeV]','[fm]','[fm]','[fm]','[fm]' do 720 i1=1,nf,1 if (incl(i1).eq.1) then write(file_out,3011) i1,lambda_f(i1),l_f(i1),j_f(i1) & ,vc_f(i1),rc_f(i1),ac_f(i1) & ,fso_f(i1)/rso_f(i1),rso_f(i1),aso_f(i1) endif 720 continue C if (case_cs.le.3) then write(file_out,*) write(file_out,3000) 'resonances in continuum states' write(file_out,3000) 'c l_f j_f E_res Gamma ' write(file_out,3000) ' [MeV] [MeV] ' do 700 i1=1,nf,1 if (incl(i1).eq.1) then do 710 i2=1,nres(i1),1 write(file_out,3001) i1,l_f(i1),j_f(i1) & ,eres(i1,i2),gres(i1,i2) 710 continue endif 700 continue endif C write(file_out,*) write(file_out,3040) '***** multipole parameters *****' write(file_out,*) write(file_out,3041) 'multipole','incl.','eff. charge' & ,'scaling factor','phase' do 730 i1=1,l_max,1 write(file_out,3042) 'E',i1,el(i1),z_eff(i1),sc(i1),phase(i1) 730 continue C write(file_out,*) write(file_out,3030) '***** scattering parameters *****' write(file_out,*) if (case_r.eq.0) then write(file_out,3077) 'non-relativistic Coulomb excitation' else write(file_out,3077) 'relativistic Coulomb excitation ' endif write(file_out,*) write(file_out,3031) 'energy of projectile',e_a,'A MeV' write(file_out,*) write(file_out,3032) 'projectile velocity',v,'c' write(file_out,3033) 'grazing angle',th_gr,'deg' write(file_out,3034) 'grazing impact parameter',b_gr,'fm' write(file_out,3035) 'grazing adiabaticity parameter',xi_min C write(file_out,*) write(file_out,3052) 'minimum scattering angle',th_min,'deg' write(file_out,3052) 'maximum scattering angle',th_max,'deg' write(file_out,3053) 'step width',dth,'deg' write(file_out,3052) 'minimum impact parameter',b_min,'fm ' C if ((case_cs.le.16).or.(case_cs.ge.19)) then write(file_out,*) write(file_out,3050) 'minimum relative energy',er_min,'MeV' write(file_out,3050) 'maximum relative energy',er_max,'MeV' write(file_out,3051) 'step width',de,'MeV' else write(file_out,*) write(file_out,3054) 'maximum relative momentum',p_max,'MeV/c' write(file_out,3055) 'step width',dp,'MeV/c' write(file_out,3050) 'maximum relative energy',er_max,'MeV' endif C write(file_out,*) C open(file_gs,status='unknown') write(file_gs,2000) '# wave function of ground state' write(file_gs,2001) '# principal quantum number n=',n_i write(file_gs,2001) '# orbital angular momentum l=',l_i write(file_gs,2002) '# spin s=',s write(file_gs,2002) '# total angular momentum j=',j_a write(file_gs,2003) '# r [fm]','f(r) [fm^(-1/2)]' do 500 i1=0,nr,1 write(file_gs,*) rv(1,i1),wf_i(n_i,i1) 500 continue close(file_gs) C open(file_sc,status='unknown') write(file_sc,2010) '# phase shifts of scattering wave functions' write(file_sc,2011) '# c',(i1,i1=1,nf) write(file_sc,2012) '# incl',(incl(i1),i1=1,nf) write(file_sc,2011) '# l',(l_f(i1),i1=1,nf) write(file_sc,2013) '# j',(j_f(i1),i1=1,nf) write(file_sc,2014) '# E',('delta',i1=1,nf) write(file_sc,2015) '# [MeV]',('[deg]',i1=1,nf) do 510 i1=0,nx,1 write(file_sc,2016) xv(i1),(delta(i2,i1),i2=1,nf) c write(30,3000) xv(i1) c & ,((dsin(delta(i2,i1)*radfac))**2,i2=1,nf) c write(31,3000) xv(i1) c & ,(ddelta(i2,i1)*radfac,i2=1,nf) 510 continue close(file_sc) C open(file_dbde,status='unknown') write(file_dbde,2020) '# reduced transition probabilities' write(file_dbde,2021) '# c',(i1,i1=1,nf) write(file_dbde,2022) '# incl',(incl(i1),i1=1,nf) write(file_dbde,2021) '# l',(l_f(i1),i1=1,nf) write(file_dbde,2023) '# j',(j_f(i1),i1=1,nf) write(file_dbde,2026) '# lambda',(lambda_f(i1),i1=1,nf) write(file_dbde,2024) '# E',('dB/dE',i1=1,nf) write(file_dbde,2025) '# [MeV]','[e^2 fm^(2*lambda)/MeV] --->' do 520 i1=0,nx,1 write(file_dbde,2027) & xv(i1),(dbde(i1,i2)/e2,i2=1,nf) 520 continue close(file_dbde) C open(file_dbdel,status='unknown') write(file_dbdel,2020) '# reduced transition probabilities' write(file_dbdel,2026) '# lambda',(i1,i1=1,l_max) write(file_dbdel,2024) '# E',('dB/dE',i1=1,l_max) write(file_dbdel,2025) '# [MeV]','[e^2 fm^(2*lambda)/MeV] --->' do 525 i1=0,nx,1 write(file_dbdel,2027) & xv(i1),(dbdel(i1,i2)/e2,i2=1,l_max) 525 continue close(file_dbdel) C open(file_ph,status='unknown') write(file_ph,2030) '# partial photo absorption cross sections' write(file_ph,2021) '# c',(i1,i1=1,nf) write(file_ph,2022) '# incl',(incl(i1),i1=1,nf) write(file_ph,2021) '# l',(l_f(i1),i1=1,nf) write(file_ph,2023) '# j',(j_f(i1),i1=1,nf) write(file_ph,2026) '# lambda',(lambda_f(i1),i1=1,nf) write(file_ph,2024) '# E',('sigma',i1=1,nf) write(file_ph,2035) '# [MeV]',('[mb]',i1=1,nf) do 530 i1=0,nx,1 write(file_ph,2027) & xv(i1),(sig_ph(i1,i2),i2=1,nf) 530 continue close(file_ph) C open(file_cp,status='unknown') write(file_cp,2040) '# partial radiative capture cross sections' write(file_cp,2021) '# c',(i1,i1=1,nf) write(file_cp,2022) '# incl',(incl(i1),i1=1,nf) write(file_cp,2021) '# l',(l_f(i1),i1=1,nf) write(file_cp,2023) '# j',(j_f(i1),i1=1,nf) write(file_cp,2026) '# lambda',(lambda_f(i1),i1=1,nf) write(file_cp,2024) '# E',('sigma',i1=1,nf) write(file_cp,2035) '# [MeV]',('[mb]',i1=1,nf) do 540 i1=0,nx,1 write(file_cp,2027) & xv(i1),(sig_cp(i1,i2),i2=1,nf) 540 continue close(file_cp) C open(file_phl,status='unknown') write(file_phl,2050) '# photo absorption cross sections' write(file_phl,2051) '#','total',(cel(i1),i1=1,l_max) write(file_phl,2052) '# E',('sigma',i1=1,(l_max+1)) write(file_phl,2053) '# [MeV]',('[mb]',i1=1,(l_max+1)) do 550 i1=0,nx,1 write(file_phl,2054) & xv(i1),sig_pht(i1),(sig_phl(i1,i2),i2=1,l_max) 550 continue close(file_phl) C open(file_cpl,status='unknown') write(file_cpl,2060) '# radiative capture cross sections' write(file_cpl,2051) '#','total',(cel(i1),i1=1,l_max) write(file_cpl,2052) '# E',('sigma',i1=1,(l_max+1)) write(file_cpl,2053) '# [MeV]',('[mb]',i1=1,(l_max+1)) do 560 i1=0,nx,1 write(file_cpl,2054) & xv(i1),sig_cpt(i1),(sig_cpl(i1,i2),i2=1,l_max) 560 continue close(file_cpl) C open(file_sf,status='unknown') if ((z_b*z_c).ne.0) then write(file_sf,2070) '# astrophysical S-factors' write(file_sf,2051) '#','total',(cel(i1),i1=1,l_max) write(file_sf,2052) '# E',(' S(E)',i1=1,(l_max+1)) write(file_sf,2073) '# [MeV]',('[MeV mb]',i1=1,(l_max+1)) do 570 i1=0,nx,1 write(file_sf,2054) & xv(i1),sfac(i1),(sfacl(i1,i2),i2=1,l_max) 570 continue else write(file_sf,*) '#no astrophysical S-factors calculated' endif close(file_sf) C if (case_cs.eq.1) then write(file_cd,2100) & '# partial Coulomb dissociation cross sections' write(file_cd,2021) '# c',(i1,i1=1,nf) write(file_cd,2022) '# incl',(incl(i1),i1=1,nf) write(file_cd,2021) '# l',(l_f(i1),i1=1,nf) write(file_cd,2023) '# j',(j_f(i1),i1=1,nf) write(file_cd,2026) '# lambda',(lambda_f(i1),i1=1,nf) write(file_cd,2101) '# E',('ds/dE',i1=1,nf) write(file_cd,2102) '# [MeV]',('[mb/MeV]',i1=1,nf) do 600 i1=0,nx,1 write(file_cd,2027) & xv(i1),(dsde(i1,i2),i2=1,nf) 600 continue write(file_cd,2111) '#',('sigma',i1=1,nf) write(file_cd,2109) '#',('mb',i1=1,nf) write(file_cd,2110) & '#',(sig(i2),i2=1,nf) C write(file_cdl,2103) '# Coulomb dissociation cross sections' write(file_cdl,2051) '#','total',(cel(i1),i1=1,l_max) write(file_cdl,2052) '# E',('ds/dE',i1=1,(l_max+1)) write(file_cdl,2104) '# [MeV]',('[mb/MeV]',i1=1,(l_max+1)) do 610 i1=0,nx,1 write(file_cdl,2054) & xv(i1),dsdet(i1),(dsdel(i1,i2),i2=1,l_max) 610 continue write(file_cdl,2112) '#',('sigma',i1=1,(l_max+1)) write(file_cdl,2113) '#',('mb',i1=1,(l_max+1)) write(file_cdl,2114) & '#',sigt,(sigl(i2),i2=1,l_max) endif C if (case_cs.eq.2) then write(file_cd,2100) & '# partial Coulomb dissociation cross sections' write(file_cd,2021) '# c',(i1,i1=1,nf) write(file_cd,2022) '# incl',(incl(i1),i1=1,nf) write(file_cd,2021) '# l',(l_f(i1),i1=1,nf) write(file_cd,2023) '# j',(j_f(i1),i1=1,nf) write(file_cd,2026) '# lambda',(lambda_f(i1),i1=1,nf) write(file_cd,2105) '# theta',('ds/dOmega',i1=1,nf) write(file_cd,2106) '# [deg]',('[mb/sr]',i1=1,nf) do 620 i1=0,nth,1 write(file_cd,2028) & thv(i1),(dsdt(i1,i2),i2=1,nf) 620 continue write(file_cd,2111) '#',('sigma',i1=1,nf) write(file_cd,2109) '#',('mb',i1=1,nf) write(file_cd,2110) & '#',(sig(i2),i2=1,nf) C write(file_cdl,2103) '# Coulomb dissociation cross sections' write(file_cdl,2051) '#','total',(cel(i1),i1=1,l_max) write(file_cdl,2107) '# theta',('ds/dOmega',i1=1,(l_max+1)) write(file_cdl,2108) '# [deg]',('[mb/sr]',i1=1,(l_max+1)) do 630 i1=0,nth,1 write(file_cdl,2055) & thv(i1),dsdtt(i1),(dsdtl(i1,i2),i2=1,l_max) 630 continue write(file_cdl,2112) '#',('sigma',i1=1,(l_max+1)) write(file_cdl,2113) '#',('mb',i1=1,(l_max+1)) write(file_cdl,2114) & '#',sigt,(sigl(i2),i2=1,l_max) endif C if (case_cs.eq.4) then write(file_cd,2100) & '# partial Coulomb dissociation cross sections' write(file_cd,2205) '# c',(i1,i1=1,nf) write(file_cd,2206) '# incl',(incl(i1),i1=1,nf) write(file_cd,2205) '# l',(l_f(i1),i1=1,nf) write(file_cd,2207) '# j',(j_f(i1),i1=1,nf) write(file_cd,2208) '# lambda',(lambda_f(i1),i1=1,nf) write(file_cd,2201) '# E',('d2s/(dEdO)',i1=1,nf) write(file_cd,2202) '# [MeV]',('[mb/(MeVsr)]',i1=1,nf) do 640 i1=0,nx,1 write(file_cd,2209) & xv(i1),(d2sde(i1,i2),i2=1,nf) 640 continue C write(file_cdl,2103) '# Coulomb dissociation cross sections' write(file_cdl,2051) '#','total',(cel(i1),i1=1,l_max) write(file_cdl,2203) '# E',('d2s/(dEdO)',i1=1,(l_max+1)) write(file_cdl,2204) '# [MeV]',('[mb/(MeVsr)]',i1=1,(l_max+1)) do 650 i1=0,nx,1 write(file_cdl,2054) & xv(i1),d2sdet(i1),(d2sdel(i1,i2),i2=1,l_max) 650 continue endif C if (case_cs.eq.5) then write(file_cd,2100) & '# partial Coulomb dissociation cross sections' write(file_cd,2205) '# c',(i1,i1=1,nf) write(file_cd,2206) '# incl',(incl(i1),i1=1,nf) write(file_cd,2205) '# l',(l_f(i1),i1=1,nf) write(file_cd,2207) '# j',(j_f(i1),i1=1,nf) write(file_cd,2208) '# lambda',(lambda_f(i1),i1=1,nf) write(file_cd,2210) '# theta',('d2s/(dEdO)',i1=1,nf) write(file_cd,2202) '# [deg]',('[mb/(MeVsr)]',i1=1,nf) do 660 i1=0,nth,1 write(file_cd,2209) & thv(i1),(d2sdt(i1,i2),i2=1,nf) 660 continue C write(file_cdl,2103) '# Coulomb dissociation cross sections' write(file_cdl,2051) '#','total',(cel(i1),i1=1,l_max) write(file_cdl,2211) '# theta',('d2s/(dEdO)',i1=1,(l_max+1)) write(file_cdl,2204) '# [deg]',('[mb/(MeVsr)]',i1=1,(l_max+1)) do 670 i1=0,nth,1 write(file_cdl,2054) & thv(i1),d2sdtt(i1),(d2sdtl(i1,i2),i2=1,l_max) 670 continue endif C if ((case_cs.eq.3).or.((case_cs.ge.6).and.(case_cs.le.8))) then write(file_cd,2300) '# Coulomb dissociation cross section' write(file_cdl,*) if (case_cs.eq.3) then write(file_cd,2301) '# theta phi ',' ds/(dO) ' write(file_cd,2301) '# [deg] [deg]',' [mb/sr] ' endif if (case_cs.eq.6) then write(file_cd,2301) '# theta phi ',' d2s/(dOdO) ' write(file_cd,2301) '# [deg] [deg]',' [mb/(sr^2)] ' endif if (case_cs.eq.7) then write(file_cd,2301) '# theta phi ',' d2s/(dEdO) ' write(file_cd,2301) '# [deg] [deg]',' [mb/(MeVsr)] ' endif if (case_cs.eq.8) then write(file_cd,2301) '# theta phi ',' d3s/(dEdOdO) ' write(file_cd,2301) '# [deg] [deg]','[mb/(MeVsr^2)]' endif do 680 i1=0,n_thbc,1 do 690 i2=0,n_phbc,1 write(file_cd,2302) thbcv(i1),phbcv(i2),sigtp(i1,i2) 690 continue 680 continue endif C if ((case_cs.ge.9).and.(case_cs.le.16)) then write(file_cd,2300) '# Coulomb dissociation cross section' write(file_cdl,*) if (case_cs.eq.9) then write(file_cd,3080) '# theta',' ds/dcos(theta) ' write(file_cd,3080) '# [deg]',' [mb] ' do 740 i1=0,n_thbc,1 write(file_cd,3081) thbcv(i1),sigtp(i1,0) 740 continue endif if (case_cs.eq.11) then write(file_cd,3080) '# theta',' ds/(dcos(theta)dO) ' write(file_cd,3080) '# [deg]',' [mb/sr] ' do 741 i1=0,n_thbc,1 write(file_cd,3081) thbcv(i1),sigtp(i1,0) 741 continue endif if (case_cs.eq.13) then write(file_cd,3080) '# theta',' d2s/(dEdcos(theta)) ' write(file_cd,3080) '# [deg]',' [mb/MeV] ' do 742 i1=0,n_thbc,1 write(file_cd,3081) thbcv(i1),sigtp(i1,0) 742 continue endif if (case_cs.eq.15) then write(file_cd,3080) '# theta','d3s/(dEdcos(theta)dO)' write(file_cd,3080) '# [deg]',' [mb/(MeV sr)] ' do 743 i1=0,n_thbc,1 write(file_cd,3081) thbcv(i1),sigtp(i1,0) 743 continue endif if (case_cs.eq.10) then write(file_cd,3080) '# phi ',' ds/dphi ' write(file_cd,3080) '# [deg]',' [mb/rad] ' do 750 i1=0,n_phbc,1 write(file_cd,3081) phbcv(i1),sigtp(0,i1) 750 continue endif if (case_cs.eq.12) then write(file_cd,3080) '# phi ',' ds/(dphidO) ' write(file_cd,3080) '# [deg]',' [mb/(rad sr)] ' do 751 i1=0,n_phbc,1 write(file_cd,3081) phbcv(i1),sigtp(0,i1) 751 continue endif if (case_cs.eq.14) then write(file_cd,3080) '# phi ',' d2s/(dEdphi) ' write(file_cd,3080) '# [deg]',' [mb/(MeV rad)] ' do 752 i1=0,n_phbc,1 write(file_cd,3081) phbcv(i1),sigtp(0,i1) 752 continue endif if (case_cs.eq.16) then write(file_cd,3080) '# phi ',' d3s/(dEdphidO) ' write(file_cd,3080) '# [deg]',' [mb/(MeV rad sr)] ' do 753 i1=0,n_phbc,1 write(file_cd,3081) phbcv(i1),sigtp(0,i1) 753 continue endif endif C if ((case_cs.eq.17).or.(case_cs.eq.18)) then write(file_cd,3060) '# longitudinal momentum distribution' write(file_cd,3061) '# p_{bc,l}','d2s/(dpldO)' if (case_cs.eq.9) & write(file_cd,3062) '# MeV/c',' [mb/((MeV/c))] ' if (case_cs.eq.10) & write(file_cd,3062) '# MeV/c','[mb/((MeV/c)sr)]' do 695 i1=(-np),np,1 write(file_cd,3063) xp(i1),dsdp(i1) 695 continue write(file_cdl,*) endif C close(file_cd) close(file_cdl) close(file_out) C 1001 format(/) 1002 format(//) 1003 format(///) 1007 format(///////) 1010 format(//////////) 1016 format(////////////////) 1028 format(////////////////////////////) C 2000 format(a31) 2001 format(a29,1x,i3) 2002 format(a29,1x,f4.1) 2003 format(a9,18x,a16) 2010 format(a43) 2011 format(a4,30(6x,i3)) 2012 format(a6,30(4x,i3,2x)) 2013 format(a4,1x,30(5x,f4.1)) 2014 format(a4,3x,30(4x,a5)) 2015 format(a7,30(4x,a5)) 2016 format(f7.4,30(1x,f8.3)) 2020 format(a34) 2021 format(a4,30(7x,i3,1x)) 2022 format(a6,30(5x,i3,3x)) 2023 format(a4,1x,30(6x,f4.1,1x)) 2024 format(a4,2x,30(5x,a5,1x)) 2025 format(a7,2x,a28) 2026 format(a8,30(3x,i3,5x)) 2027 format(f7.4,30(1x,e10.4)) 2028 format(f7.3,30(1x,e10.4)) 2030 format(a41) 2035 format(a7,30(5x,a4,2x)) 2040 format(a42) 2050 format(a33) 2051 format(a1,10x,a5,30(9x,a3,1x)) 2052 format(a4,3x,31(4x,a5,4x)) 2053 format(a7,31(5x,a4,4x)) 2054 format(f7.4,31(1x,e12.6)) 2055 format(f7.3,31(1x,e12.6)) 2060 format(a34) 2070 format(a25) 2073 format(a7,31(3x,a8,2x)) 2100 format(a45) 2101 format(a4,2x,30(5x,a5,1x)) 2102 format(a7,30(3x,a8)) 2103 format(a36) 2104 format(a7,31(3x,a8,2x)) 2105 format(a7,30(1x,a9,1x)) 2106 format(a7,30(3x,a7,1x)) 2107 format(a7,31(2x,a9,2x)) 2108 format(a7,30(3x,a7,3x)) 2109 format(a,6x,30(6x,a2,3x)) 2110 format(a1,6x,30(1x,e10.4)) 2111 format(a1,5x,30(5x,a5,1x)) 2112 format(a1,6x,31(4x,a5,4x)) 2113 format(a1,6x,31(6x,a2,5x)) 2114 format(a1,6x,31(1x,e12.6)) 2201 format(a4,2x,30(3x,a10)) 2202 format(a7,30(1x,a12)) 2203 format(a4,4x,31(1x,a10,2x)) 2204 format(a7,31(1x,a12)) 2205 format(a4,30(8x,i3,2x)) 2206 format(a6,30(6x,i3,4x)) 2207 format(a4,1x,30(7x,f4.1,2x)) 2208 format(a8,30(4x,i3,6x)) 2209 format(f7.4,30(2x,e11.5)) 2210 format(a7,30(2x,a10,1x)) 2211 format(a7,1x,31(1x,a10,2x)) 2300 format(a36) 2301 format(a17,3x,a14) 2302 format(1x,f7.3,2x,f7.3,4x,e11.5) C 3000 format(4x,a30) 3001 format(3x,i2,2x,i2,1x,f5.1,1x,f7.4,1x,f8.4) 3002 format(3x,a38) 3003 format(4x,a3,i2,2x,a3,i2,2x,a3,f4.1,2x,a3,1x,f4.1) 3004 format(4x,a14,1x,f8.4,1x,a3) 3005 format(4x,a20) 3006 format(7x,a3,3x,a3,3x,a3,4x,a4,3x,a4,2x,a4) 3007 format(6x,a5,2x,a4,2x,a4,10x,a4,2x,a4) 3008 format(4x,f7.4,1x,f5.3,1x,f5.3,1x,f7.5,1x,f5.3,1x,f5.3) 3009 format(3x,a28) 3010 format(4x,a,1x,a6,1x,a,3x,a,5x,a3,4x,a3,3x,a3,4x,a4,3x,a4,2x,a4) 3011 format(3x,i2,4x,i1,2x,i2,1x,f4.1,2x,f7.4,1x,f5.3,1x,f5.3 & ,1x,f7.5,1x,f5.3,1x,f5.3) 3012 format(23x,a5,2x,a4,2x,a4,10x,a4,2x,a4) 3013 format(3x,a32) 3014 format(4x,a6,8x,a3,i3,3x,a3,i4) 3015 format(4x,f10.6,1x,a1,3x,f11.3,1x,a3) 3016 format(4x,a10,4x,a3,i3,2x,a3,i4) 3020 format(3x,a32) 3030 format(3x,a33) 3031 format(4x,a20,10x,f8.3,1x,a5) 3032 format(4x,a19,14x,f7.5,1x,a1) 3033 format(4x,a13,19x,f6.3,1x,a3) 3034 format(4x,a24,7x,f7.3,1x,a2) 3035 format(4x,a30,2x,f6.3) 3040 format(3x,a32) 3041 format(4x,a9,3x,a5,3x,a11,3x,a14,3x,a5) 3042 format(7x,a1,i2,8x,i1,6x,f8.5,8x,f7.3,6x,f7.3) 3050 format(4x,a23,4x,f7.3,1x,a3) 3051 format(4x,a10,17x,f7.3,1x,a3) 3052 format(4x,a24,3x,f7.3,1x,a3) 3053 format(4x,a10,17x,f7.3,1x,a3) 3054 format(4x,a25,2x,f7.3,1x,a5) 3055 format(4x,a10,17x,f7.3,1x,a5) 3060 format(a36) 3061 format(a10,8x,a11) 3062 format(a8,7x,a16) 3063 format(1x,f8.3,7x,e13.6) 3070 format(4x,a12) 3071 format(4x,a3,i3,5x,a3,i4,5x,a3,f4.1) 3072 format(4x,a3,f9.5,1x,a4,f10.3,1x,a3) 3073 format(4x,a21) 3074 format(4x,a28,1x,f5.2,1x,a2) 3075 format(4x,a14,13x,f7.2,1x,a2) 3076 format(4x,a16,11x,i4) 3077 format(4x,a35) 3080 format(a7,2x,a21) 3081 format(1x,f7.3,6x,e11.5) C stop end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ***************************************************************** subroutine init_simul(e_i,de,dth,z_x,case_r,nf,nx,nth,file,lamax) C ***************************************************************** implicit none integer*4 case,case_r,nth_max,nx_max,l_max,tl_max,n_max parameter (nth_max=90,nx_max=501,n_max=5000, & l_max=3,tl_max=2*l_max) integer*4 i1,i2,l,m,nf,nx,nth,z_x,n_thbc,n_phbc,file,lamax real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,e_i,e_x,a,v,gamma,p_bc,theta,th_bc,ph_bc,dthr & ,tmp0,tmp1,tmp2,mu_bc,h2m,zze2,d3s,res,de,dth & ,sigtp(0:nth_max,0:nth_max) & ,xv(0:nx_max),wf_f(0:n_max) & ,thv(0:nth_max) complex*16 i,wffac & ,coeff_a(0:tl_max,-tl_max:tl_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /d3s_cd/ sigtp common /wave_f/ xv,wf_f,wffac common /par_bc/ mu_bc,h2m,zze2 common /angles/ thv common /co_a/ coeff_a C tmp0 = fpi*dfloat(z_x)/(2.d00*hbarc) tmp0 = tmp0*tmp0*e2*mu_bc tmp1 = tpi*hbarc tmp0 = tmp0/(tmp1*tmp1*tmp1) C factor 10: fm^2 -> mb tmp0 = tmp0*10.d00 C do 100 i1=0,nth,1 theta = thv(i1) do 110 i2=0,nx,1 e_x = e_i+xv(i2) p_bc = dsqrt(2.d00*mu_bc*xv(i2)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ C ++++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,i2,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif write(file,*) tmp2 do 120 l=0,tl_max,1 do 130 m=(-l),l,1 write(file,*) coeff_a(l,m) 130 continue 120 continue 110 continue 100 continue C return end C ************************* subroutine get_say(x,say) C ************************* implicit none integer*4 l_max,tl_max,i1,i2 parameter (l_max=3,tl_max=2*l_max) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,x,say,plx(0:tl_max) complex*16 i,ctmp & ,coeff_a(0:tl_max,-tl_max:tl_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /co_a/ coeff_a C plx(0) = 1.d00 plx(1) = x do 120 i1=2,tl_max,1 plx(i1) = (dfloat(2*i1-1)*x*plx(i1-1) & -dfloat(i1-1)*plx(i1-2))/dfloat(i1) 120 continue C ctmp = 0.d00 do 100 i1=0,tl_max,1 ctmp = ctmp+coeff_a(i1,0)*plx(i1)*dsqrt(dfloat(2*i1+1)) 100 continue C say = dreal(ctmp)/dsqrt(fpi) C return end C ************************************************* subroutine get_dsdp(e_i,de,dp,dth & ,z_x,case_r,nf,nx,np,nth,case,lamax) C ************************************************* implicit none integer*4 nx_max,n_max,nth_max parameter (nx_max=501,n_max=5000,nth_max=90) integer*4 case,case_r,nf,nx,np,nth,z_x,i1,i2,i3,nn,lamax real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,mu_bc,h2m,zze2,e_i,de,dp,dth,theta,e_x,dthr & ,a,v,gamma,x,say,res,dsdp(-nx_max:nx_max) & ,tmp0,tmp1,tmp2,tmp3(0:n_max),tmp4(0:nth_max) & ,tmp5(0:nth_max,-nx_max:nx_max),xp(-nx_max:nx_max) & ,xv(0:nx_max),wf_f(0:n_max),thv(0:nth_max) complex*16 i,wffac C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /wave_f/ xv,wf_f,wffac common /par_bc/ mu_bc,h2m,zze2 common /angles/ thv common /momenta/ xp common /momdistr/ dsdp C tmp0 = fpi*dfloat(z_x)/(2.d00*hbarc) tmp0 = tmp0*tmp0*e2 tmp1 = tpi*hbarc tmp0 = tpi*tmp0/(tmp1*tmp1*tmp1) C factor 10: fm^2 -> mb tmp0 = 10.d00*tmp0 C if (case.eq.1) then dthr = dth*radfac do 200 i3=0,nth,1 tmp4(i3) = dsin(radfac*thv(i3)) 200 continue tmp0 = tmp0*tpi do 210 i3=0,nth,1 theta = thv(i3) if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp0/(tmp1*tmp1*tmp1*tmp1) do 220 i1=(-np),np,1 nn = np-iabs(i1) do 230 i2=iabs(i1),np,1 e_x = e_i+xv(i2) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ C ++++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,i2,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++ tmp2 = a*a*xp(i2) if (i2.ne.0) then x = xp(i1)/xp(iabs(i2)) else x = 0.d00 endif C +++++++++++++++++++ call get_say(x,say) C +++++++++++++++++++ tmp3(i2-iabs(i1)) = say*tmp2 230 continue C ++++++++++++++++++++++++++++ call simpson(nn,tmp3,dp,res) C ++++++++++++++++++++++++++++ tmp5(i3,i1) = res*tmp1 220 continue else do 240 i1=(-np),np,1 tmp5(i3,i1) = 0.d00 240 continue endif 210 continue do 250 i1=(-np),np,1 do 260 i3=0,nth,1 tmp3(i3) = tmp5(i3,i1)*tmp4(i3) 260 continue C +++++++++++++++++++++++++++++++ call simpson(nth,tmp3,dthr,res) C +++++++++++++++++++++++++++++++ dsdp(i1) = res 250 continue endif C if (case.eq.2) then theta = thv(0) if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp0/(tmp1*tmp1*tmp1*tmp1) do 100 i1=(-np),np,1 nn = np-iabs(i1) do 110 i2=iabs(i1),np,1 e_x = e_i+xv(i2) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ C ++++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,i2,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++ tmp2 = a*a*xp(i2) if (i2.ne.0) then x = xp(i1)/xp(iabs(i2)) else x = 0.d00 endif C +++++++++++++++++++ call get_say(x,say) C +++++++++++++++++++ tmp3(i2-iabs(i1)) = say*tmp2 110 continue C ++++++++++++++++++++++++++++ call simpson(nn,tmp3,dp,res) C ++++++++++++++++++++++++++++ dsdp(i1) = res*tmp1 100 continue else do 120 i1=(-np),np,1 dsdp(i1) = 0.d00 120 continue endif endif C return end C *********************************** subroutine get_resonances(de,nf,nx) C *********************************** implicit none integer*4 nx_max,n_max,nf_max,nres_max,l_max parameter (nx_max=501,n_max=5000,nf_max=20,nres_max=5,l_max=3) integer*4 nf,nx,i1,i2,nres(nf_max) & ,l_f(nf_max),incl(nf_max),el(0:l_max) & ,lambda_f(nf_max) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,de,tmp1,tmp2,tmp3,tmp4,aa,bb,cc & ,eres(nf_max,nres_max),gres(nf_max,nres_max) & ,xv(0:nx_max),wf_f(0:n_max),j_f(nf_max) & ,delta(nf_max,0:nx_max),ddelta(nf_max,0:nx_max) complex*16 i,wffac,z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /phases/ delta,ddelta common /wave_f/ xv,wf_f,wffac common /reson/ eres,gres,nres common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el C tmp1 = 2.d00*de tmp2 = tmp1*de do 100 i1=1,nf,1 if (incl(i1).eq.1) then do 110 i2=1,(nx-1),1 ddelta(i1,i2) = (delta(i1,i2+1)-delta(i1,i2-1))/tmp1 c write(40,*) i1,i2,delta(i1,i2),ddelta(i1,i2) 110 continue endif 100 continue C do 120 i1=1,nf,1 nres(i1) = 0 if (incl(i1).eq.1) then do 130 i2=2,nx-2,1 if ((ddelta(i1,i2).gt.ddelta(i1,i2-1)) & .and.(ddelta(i1,i2).gt.ddelta(i1,i2+1))) then if (nres(i1).lt.nres_max) then nres(i1) = nres(i1)+1 aa = ddelta(i1,i2) bb = (ddelta(i1,i2+1)-ddelta(i1,i2-1))/tmp1 cc = (ddelta(i1,i2+1)+ddelta(i1,i2-1) & -2.d00*ddelta(i1,i2))/tmp2 tmp3 = -bb/(2.d00*cc) eres(i1,nres(i1)) = xv(i2)+tmp3 tmp4 = aa+(bb+cc*tmp3)*tmp3 gres(i1,nres(i1)) = 2.d00/(radfac*tmp4) endif endif 130 continue endif 120 continue C return end C ************************************************************* subroutine get_dxscd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,case,case_tp,lamax) C ************************************************************* implicit none integer*4 case,case_r,case_tp,nth_max,nx_max,n_max parameter (nth_max=90,nx_max=501,n_max=5000) integer*4 i1,i2,i3,i4,nf,nx,nth,z_x,n_thbc,n_phbc,lamax real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,e_i,e_x,a,v,gamma,p_bc,theta,th_bc,ph_bc,dthr & ,tmp0,tmp1,tmp2,mu_bc,h2m,zze2,d3s,res,de,dth & ,sigtp(0:nth_max,0:nth_max),tmp5(0:nth_max) & ,xv(0:nx_max),wf_f(0:n_max),tmp6(0:nx_max) & ,thbcv(0:nth_max),phbcv(0:nth_max),thv(0:nth_max) & ,tmp3(0:nth_max,0:nx_max),tmp4(0:n_max) complex*16 i,wffac C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /d3s_cd/ sigtp common /wave_f/ xv,wf_f,wffac common /par_bc/ mu_bc,h2m,zze2 common /angles/ thv common /angbc/ thbcv,phbcv C tmp0 = fpi*dfloat(z_x)/(2.d00*hbarc) tmp0 = tmp0*tmp0*e2*mu_bc tmp1 = tpi*hbarc tmp0 = tmp0/(tmp1*tmp1*tmp1) C factor 10: fm^2 -> mb tmp0 = tmp0*10.d00 C dthr = dth*radfac C if (case.eq.0) then do 480 i4=0,nth,1 tmp5(i4) = dsin(radfac*thv(i4)) 480 continue do 440 i1=0,n_thbc,1 th_bc = thbcv(i1) do 420 i2=0,n_phbc,1 ph_bc = phbcv(i2) do 400 i3=0,nx,1 e_x = e_i+xv(i3) p_bc = dsqrt(2.d00*mu_bc*xv(i3)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ do 470 i4=0,nth,1 theta = thv(i4) C ++++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,i3,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif C +++++++++++++++++++++++++++++++++++++ call get_dxs(th_bc,ph_bc,d3s,case_tp) C +++++++++++++++++++++++++++++++++++++ tmp4(i4) = d3s*tmp2*tmp5(i4) 470 continue C +++++++++++++++++++++++++++++++ call simpson(nth,tmp4,dthr,res) C +++++++++++++++++++++++++++++++ tmp6(i3) = res 400 continue do 490 i3=0,nx,1 tmp4(i3) = tmp6(i3) 490 continue C ++++++++++++++++++++++++++++ call simpson(nx,tmp4,de,res) C ++++++++++++++++++++++++++++ sigtp(i1,i2) = res*tpi 420 continue 440 continue endif C if (case.eq.1) then theta = thv(0) do 140 i1=0,n_thbc,1 th_bc = thbcv(i1) do 100 i3=0,nx,1 e_x = e_i+xv(i3) p_bc = dsqrt(2.d00*mu_bc*xv(i3)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ C ++++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,i3,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif do 120 i2=0,n_phbc,1 ph_bc = phbcv(i2) C +++++++++++++++++++++++++++++++++++++ call get_dxs(th_bc,ph_bc,d3s,case_tp) C +++++++++++++++++++++++++++++++++++++ tmp3(i2,i3) = d3s*tmp2 120 continue 110 continue 100 continue do 150 i2=0,n_phbc,1 do 160 i3=0,nx,1 tmp4(i3) = tmp3(i2,i3) 160 continue C ++++++++++++++++++++++++++++ call simpson(nx,tmp4,de,res) C ++++++++++++++++++++++++++++ sigtp(i1,i2) = res 150 continue 140 continue endif C if (case.eq.2) then e_x = e_i+xv(0) p_bc = dsqrt(2.d00*mu_bc*xv(0)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ do 270 i3=0,nth,1 tmp5(i3) = dsin(radfac*thv(i3)) 270 continue do 240 i1=0,n_thbc,1 th_bc = thbcv(i1) do 200 i3=0,nth,1 theta = thv(i3) C +++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,0,lamax) C +++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif do 220 i2=0,n_phbc,1 ph_bc = phbcv(i2) C +++++++++++++++++++++++++++++++++++++ call get_dxs(th_bc,ph_bc,d3s,case_tp) C +++++++++++++++++++++++++++++++++++++ tmp3(i2,i3) = d3s*tmp2 220 continue 210 continue 200 continue do 250 i2=0,n_phbc,1 do 260 i3=0,nth,1 tmp4(i3) = tmp3(i2,i3)*tmp5(i3) 260 continue C +++++++++++++++++++++++++++++++ call simpson(nth,tmp4,dthr,res) C +++++++++++++++++++++++++++++++ sigtp(i1,i2) = res*tpi 250 continue 240 continue endif C if (case.eq.3) then e_x = e_i+xv(0) p_bc = dsqrt(2.d00*mu_bc*xv(0)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ theta = thv(0) C +++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,0,lamax) C +++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif C do 310 i1=0,n_thbc,1 th_bc = thbcv(i1) do 320 i2=0,n_phbc,1 ph_bc = phbcv(i2) C +++++++++++++++++++++++++++++++++++++ call get_dxs(th_bc,ph_bc,d3s,case_tp) C +++++++++++++++++++++++++++++++++++++ sigtp(i1,i2) = d3s*tmp2 320 continue 310 continue endif C return end C ************************************************************* subroutine get_d3scd(e_i,de,dth & ,z_x,case_r,nf,nx,nth,n_thbc,n_phbc,case,lamax) C ************************************************************* implicit none integer*4 case,case_r,nth_max,nx_max,n_max parameter (nth_max=90,nx_max=501,n_max=5000) integer*4 i1,i2,i3,i4,nf,nx,nth,z_x,n_thbc,n_phbc,lamax real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,e_i,e_x,a,v,gamma,p_bc,theta,th_bc,ph_bc,dthr & ,tmp0,tmp1,tmp2,mu_bc,h2m,zze2,d3s,res,de,dth & ,sigtp(0:nth_max,0:nth_max),tmp5(0:nth_max) & ,xv(0:nx_max),wf_f(0:n_max),tmp6(0:nx_max) & ,thbcv(0:nth_max),phbcv(0:nth_max),thv(0:nth_max) & ,tmp3(0:nth_max,0:nx_max),tmp4(0:n_max) complex*16 i,wffac C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /d3s_cd/ sigtp common /wave_f/ xv,wf_f,wffac common /par_bc/ mu_bc,h2m,zze2 common /angles/ thv common /angbc/ thbcv,phbcv C tmp0 = fpi*dfloat(z_x)/(2.d00*hbarc) tmp0 = tmp0*tmp0*e2*mu_bc tmp1 = tpi*hbarc tmp0 = tmp0/(tmp1*tmp1*tmp1) C factor 10: fm^2 -> mb tmp0 = tmp0*10.d00 C dthr = dth*radfac C if (case.eq.0) then do 480 i4=0,nth,1 tmp5(i4) = dsin(radfac*thv(i4)) 480 continue do 440 i1=0,n_thbc,1 th_bc = thbcv(i1) do 420 i2=0,n_phbc,1 ph_bc = phbcv(i2) do 400 i3=0,nx,1 e_x = e_i+xv(i3) p_bc = dsqrt(2.d00*mu_bc*xv(i3)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ do 470 i4=0,nth,1 theta = thv(i4) C ++++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,i3,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif C +++++++++++++++++++++++++++++ call get_d3s(th_bc,ph_bc,d3s) C +++++++++++++++++++++++++++++ tmp4(i4) = d3s*tmp2*tmp5(i4) 470 continue C +++++++++++++++++++++++++++++++ call simpson(nth,tmp4,dthr,res) C +++++++++++++++++++++++++++++++ tmp6(i3) = res 400 continue do 490 i3=0,nx,1 tmp4(i3) = tmp6(i3) 490 continue C ++++++++++++++++++++++++++++ call simpson(nx,tmp4,de,res) C ++++++++++++++++++++++++++++ sigtp(i1,i2) = res*tpi 420 continue 440 continue endif C if (case.eq.1) then theta = thv(0) do 140 i1=0,n_thbc,1 th_bc = thbcv(i1) do 100 i3=0,nx,1 e_x = e_i+xv(i3) p_bc = dsqrt(2.d00*mu_bc*xv(i3)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ C ++++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,i3,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif do 120 i2=0,n_phbc,1 ph_bc = phbcv(i2) C +++++++++++++++++++++++++++++ call get_d3s(th_bc,ph_bc,d3s) C +++++++++++++++++++++++++++++ tmp3(i2,i3) = d3s*tmp2 120 continue 110 continue 100 continue do 150 i2=0,n_phbc,1 do 160 i3=0,nx,1 tmp4(i3) = tmp3(i2,i3) 160 continue C ++++++++++++++++++++++++++++ call simpson(nx,tmp4,de,res) C ++++++++++++++++++++++++++++ sigtp(i1,i2) = res 150 continue 140 continue endif C if (case.eq.2) then e_x = e_i+xv(0) p_bc = dsqrt(2.d00*mu_bc*xv(0)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ do 270 i3=0,nth,1 tmp5(i3) = dsin(radfac*thv(i3)) 270 continue do 240 i1=0,n_thbc,1 th_bc = thbcv(i1) do 200 i3=0,nth,1 theta = thv(i3) C +++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,0,lamax) C +++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif do 220 i2=0,n_phbc,1 ph_bc = phbcv(i2) C +++++++++++++++++++++++++++++ call get_d3s(th_bc,ph_bc,d3s) C +++++++++++++++++++++++++++++ tmp3(i2,i3) = d3s*tmp2 220 continue 210 continue 200 continue do 250 i2=0,n_phbc,1 do 260 i3=0,nth,1 tmp4(i3) = tmp3(i2,i3)*tmp5(i3) 260 continue C +++++++++++++++++++++++++++++++ call simpson(nth,tmp4,dthr,res) C +++++++++++++++++++++++++++++++ sigtp(i1,i2) = res*tpi 250 continue 240 continue endif C if (case.eq.3) then e_x = e_i+xv(0) p_bc = dsqrt(2.d00*mu_bc*xv(0)) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ theta = thv(0) C +++++++++++++++++++++++++++++++++++++++++++++++++ call get_coeff_a(e_x,v,a,theta,case_r,nf,0,lamax) C +++++++++++++++++++++++++++++++++++++++++++++++++ if (theta.gt.0.d00) then tmp1 = dsin(theta*radfac/2.d00) tmp1 = tmp1*tmp1 tmp2 = tmp0*a*a*p_bc/(tmp1*tmp1) else tmp2 = 0.d00 endif C do 310 i1=0,n_thbc,1 th_bc = thbcv(i1) do 320 i2=0,n_phbc,1 ph_bc = phbcv(i2) C +++++++++++++++++++++++++++++ call get_d3s(th_bc,ph_bc,d3s) C +++++++++++++++++++++++++++++ sigtp(i1,i2) = d3s*tmp2 320 continue 310 continue endif C return end C **************************************************** subroutine get_coeff_a(e,v,a,theta,case,nf,ne,lamax) C **************************************************** implicit none integer*4 nf_max,l_max,nx_max,tl_max parameter (nf_max=20,l_max=3,tl_max=2*l_max,nx_max=501) integer*4 case,nf,ne,i1,i2,i3,i4,i5,i6,lamax & ,l_f(nf_max),incl(nf_max),el(0:l_max) & ,lambda_f(nf_max) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,e,v,a,theta,j_f(nf_max) & ,coeff_b(nf_max,nf_max,0:(2*l_max)) complex*16 i,ctmp,rmat(0:nx_max,nf_max) & ,telm(l_max,l_max,0:tl_max,-tl_max:tl_max) & ,coeff_a(0:tl_max,-tl_max:tl_max),z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /co_b/ coeff_b common /redmat/ rmat common /excf2/ telm common /co_a/ coeff_a C C +++++++++++++++++++++++++++++++++++++ call get_telm(e,v,a,theta,case,lamax) C +++++++++++++++++++++++++++++++++++++ C do 300 i3=0,tl_max,1 do 310 i4=(-tl_max),tl_max,1 coeff_a(i3,i4) = 0.d00 310 continue 300 continue C do 100 i1=1,nf,1 if (incl(i1).eq.1) then i5 = lambda_f(i1) if (el(i5).eq.1) then do 110 i2=1,nf,1 if (incl(i2).eq.1) then i6 = lambda_f(i2) if (el(i6).eq.1) then do 200 i3=0,tl_max,1 ctmp = coeff_b(i1,i2,i3) & *rmat(ne,i1)*dconjg(rmat(ne,i2)) do 210 i4=(-i3),i3,1 coeff_a(i3,i4) = coeff_a(i3,i4)+ctmp*telm(i5,i6,i3,i4) 210 continue 200 continue endif endif 110 continue endif endif 100 continue C c c do 400 i3=0,tl_max,1 c do 410 i4=(-i3),i3,1 c write(40,*) i3,i4,coeff_a(i3,i4) c 410 continue c 400 continue C return end C ***************************************** subroutine get_dxs(theta,phi,d3s,case_tp) C ***************************************** implicit none integer*4 l_max,tl_max,i1,i2,case_tp parameter (l_max=3,tl_max=2*l_max) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,theta,phi,d3s complex*16 i,ctmp,ylm(0:tl_max,-tl_max:tl_max) & ,coeff_a(0:tl_max,-tl_max:tl_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /co_a/ coeff_a common /ylmtp/ ylm C C ++++++++++++++++++++++++++++++++ call get_ylmx(theta,phi,case_tp) C ++++++++++++++++++++++++++++++++ C ctmp = 0.d00 do 100 i1=0,tl_max,1 do 110 i2=(-i1),i1,1 ctmp = ctmp+coeff_a(i1,i2)*ylm(i1,i2) 110 continue 100 continue C d3s = dreal(ctmp) C return end C ********************************* subroutine get_d3s(theta,phi,d3s) C ********************************* implicit none integer*4 l_max,tl_max,i1,i2 parameter (l_max=3,tl_max=2*l_max) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,theta,phi,d3s complex*16 i,ctmp,ylm(0:tl_max,-tl_max:tl_max) & ,coeff_a(0:tl_max,-tl_max:tl_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /co_a/ coeff_a common /ylmtp/ ylm C C +++++++++++++++++++++++ call get_ylm(theta,phi) C +++++++++++++++++++++++ C ctmp = 0.d00 do 100 i1=0,tl_max,1 do 110 i2=(-i1),i1,1 ctmp = ctmp+coeff_a(i1,i2)*ylm(i1,i2) 110 continue 100 continue C d3s = dreal(ctmp) C return end C ****************************************************** subroutine get_d2sdet(e_i,de,dth & ,nf,nx,nth,z_x,case_r,case,lamax) C ****************************************************** implicit none integer*4 nf,nx,nth,z_x,case_r,case & ,nth_max,nx_max,n_max,nf_max,l_max parameter (nth_max=90,nx_max=501,n_max=5000,nf_max=20 & ,l_max=3) integer*4 l_f(nf_max),incl(nf_max),el(0:l_max) & ,lambda_f(nf_max),i1,i2,i3,lamax real*8 e_i,de,dth,dthr,e_x,a,v,gamma,theta,thv(0:nth_max) & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,xv(0:nx_max),wf_f(0:n_max) & ,tmp0(0:nth_max),tmp1(0:n_max) & ,tmp2(0:nx_max,0:nth_max,nf_max) & ,res,j_f(nf_max) & ,d2sde(0:nx_max,nf_max) & ,d2sdel(0:nx_max,0:l_max),d2sdet(0:nx_max) & ,d2sdt(0:nth_max,nf_max) & ,d2sdtl(0:nth_max,0:l_max),d2sdtt(0:nth_max) & ,sig(nf_max),sigl(0:l_max),sigt & ,d2s(nf_max) complex*16 i,wffac,z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /wave_f/ xv,wf_f,wffac common /angles/ thv common /sig_cd/ sig,sigl,sigt common /d2s_cd/ d2s common /d2s_cde/ d2sde,d2sdel,d2sdet common /d2s_cdt/ d2sdt,d2sdtl,d2sdtt C if (case.eq.1) then theta = thv(0) do 100 i1=0,nx,1 e_x = e_i+xv(i1) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ C +++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d2s_cd(e_x,a,v,theta,i1,z_x,case_r,nf,lamax) C +++++++++++++++++++++++++++++++++++++++++++++++++++++ do 110 i2=1,nf,1 if (incl(i2).eq.1) then i3 = lambda_f(i2) d2sde(i1,i2) = d2s(i2) else d2sde(i1,i2) = 0.d00 endif 110 continue do 120 i2=0,l_max,1 d2sdel(i1,i2) = 0.d00 120 continue d2sdet(i1) = 0.d00 do 130 i2=1,nf,1 i3 = lambda_f(i2) if (el(i3).eq.1) then d2sdel(i1,i3) = d2sdel(i1,i3)+d2sde(i1,i2) d2sdet(i1) = d2sdet(i1)+d2sde(i1,i2) endif 130 continue 100 continue endif C if (case.eq.2) then e_x = e_i+xv(0) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ do 200 i1=0,nth,1 theta = thv(i1) C ++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d2s_cd(e_x,a,v,theta,0,z_x,case_r,nf,lamax) C ++++++++++++++++++++++++++++++++++++++++++++++++++++ do 210 i2=1,nf,1 if (incl(i2).eq.1) then i3 = lambda_f(i2) d2sdt(i1,i2) = d2s(i2) else d2sdt(i1,i2) = 0.d00 endif 210 continue do 220 i2=0,l_max,1 d2sdtl(i1,i2) = 0.d00 220 continue d2sdtt(i1) = 0.d00 do 230 i2=1,nf,1 i3 = lambda_f(i2) if (el(i3).eq.1) then d2sdtl(i1,i3) = d2sdtl(i1,i3)+d2sdt(i1,i2) d2sdtt(i1) = d2sdtt(i1)+d2sdt(i1,i2) endif 230 continue 200 continue endif C return end C **************************************************************** subroutine get_dsdet(e_i,de,dth,nf,nx,nth,z_x,case_r,case,lamax) C **************************************************************** implicit none integer*4 nf,nx,nth,z_x,case_r,case & ,nth_max,nx_max,n_max,nf_max,l_max parameter (nth_max=90,nx_max=501,n_max=5000,nf_max=20 & ,l_max=3) integer*4 l_f(nf_max),incl(nf_max),el(0:l_max) & ,lambda_f(nf_max),i1,i2,i3,lamax real*8 e_i,de,dth,dthr,e_x,a,v,gamma,theta,thv(0:nth_max) & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,xv(0:nx_max),wf_f(0:n_max) & ,tmp0(0:nth_max),tmp1(0:n_max) & ,tmp2(0:nx_max,0:nth_max,nf_max) & ,res,j_f(nf_max) & ,dsde(0:nx_max,nf_max) & ,dsdel(0:nx_max,0:l_max),dsdet(0:nx_max) & ,dsdt(0:nth_max,nf_max) & ,dsdtl(0:nth_max,0:l_max),dsdtt(0:nth_max) & ,sig(nf_max),sigl(0:l_max),sigt & ,d2s(nf_max) complex*16 i,wffac,z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /wave_f/ xv,wf_f,wffac common /angles/ thv common /sig_cd/ sig,sigl,sigt common /d2s_cd/ d2s common /ds_cde/ dsde,dsdel,dsdet common /ds_cdt/ dsdt,dsdtl,dsdtt C do 100 i1=0,nx,1 e_x = e_i+xv(i1) C ++++++++++++++++++++++++++++++ call kinematics(e_x,a,v,gamma) C ++++++++++++++++++++++++++++++ do 110 i2=0,nth,1 theta = thv(i2) C +++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d2s_cd(e_x,a,v,theta,i1,z_x,case_r,nf,lamax) C +++++++++++++++++++++++++++++++++++++++++++++++++++++ do 120 i3=1,nf,1 if (incl(i3).eq.1) then tmp2(i1,i2,i3) = d2s(i3) else tmp2(i1,i2,i3) = 0.d00 endif 120 continue 110 continue 100 continue C dthr = dth*radfac do 200 i1=0,nth,1 tmp0(i1) = dsin(thv(i1)*radfac) 200 continue C sigt = 0.d00 do 420 i2=0,l_max,1 sigl(i2) = 0.d00 420 continue C if (case.eq.1) then do 240 i3=1,nf,1 do 250 i1=0,nx,1 dsde(i1,i3) = 0.d00 250 continue 240 continue do 210 i3=1,nf,1 if (incl(i3).eq.1) then do 220 i1=0,nx,1 do 230 i2=0,nth,1 tmp1(i2) = tmp2(i1,i2,i3)*tmp0(i2) 230 continue C +++++++++++++++++++++++++++++++ call simpson(nth,tmp1,dthr,res) C +++++++++++++++++++++++++++++++ dsde(i1,i3) = tpi*res 220 continue endif 210 continue do 260 i1=0,nx,1 do 280 i2=0,l_max,1 dsdel(i1,i2) = 0.d00 280 continue dsdet(i1) = 0.d00 do 270 i3=1,nf,1 i2 = lambda_f(i3) if (el(i2).eq.1) then dsdel(i1,i2) = dsdel(i1,i2)+dsde(i1,i3) dsdet(i1) = dsdet(i1)+dsde(i1,i3) endif 270 continue 260 continue C do 400 i3=1,nf,1 if (incl(i3).eq.1) then i2 = lambda_f(i3) do 410 i1=0,nx,1 tmp1(i1) = dsde(i1,i3) 410 continue C ++++++++++++++++++++++++++++ call simpson(nx,tmp1,de,res) C ++++++++++++++++++++++++++++ sig(i3) = res sigl(i2) = sigl(i2)+sig(i3) sigt = sigt+sig(i3) else sig(i3) = 0.d00 endif 400 continue C endif C if (case.eq.2) then do 340 i3=1,nf,1 do 350 i1=0,nth,1 dsdt(i1,i3) = 0.d00 350 continue 340 continue do 310 i3=1,nf,1 if (incl(i3).eq.1) then do 320 i1=0,nth,1 do 330 i2=0,nx,1 tmp1(i2) = tmp2(i2,i1,i3) 330 continue C ++++++++++++++++++++++++++++ call simpson(nx,tmp1,de,res) C ++++++++++++++++++++++++++++ dsdt(i1,i3) = res 320 continue endif 310 continue do 360 i1=0,nth,1 do 380 i2=0,l_max,1 dsdtl(i1,i2) = 0.d00 380 continue dsdtt(i1) = 0.d00 do 370 i3=1,nf,1 i2 = lambda_f(i3) if (el(i2).eq.1) then dsdtl(i1,i2) = dsdtl(i1,i2)+dsdt(i1,i3) dsdtt(i1) = dsdtt(i1)+dsdt(i1,i3) endif 370 continue 360 continue C do 500 i3=1,nf,1 if (incl(i3).eq.1) then i2 = lambda_f(i3) do 510 i1=0,nth,1 tmp1(i1) = dsdt(i1,i3)*tmp0(i1) 510 continue C +++++++++++++++++++++++++++++++ call simpson(nth,tmp1,dthr,res) C +++++++++++++++++++++++++++++++ sig(i3) = tpi*res sigl(i2) = sigl(i2)+sig(i3) sigt = sigt+sig(i3) else sig(i3) = 0.d00 endif 500 continue C endif C return end C ***************************************************** subroutine get_d2s_cd(e,a,v,theta,ne,z,case,nf,lamax) C ***************************************************** implicit none integer*4 case,l_max,tl_max,nf_max,z,ne,nx_max,nf parameter (l_max=3,tl_max=2*l_max,nf_max=20,nx_max=501) integer*4 l_f(nf_max),incl(nf_max),el(0:l_max) & ,lambda_f(nf_max),i1,i2,lamax real*8 e,a,v,theta,tmp1,tmp2(0:l_max) & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,dbde(0:nx_max,nf_max),dbdel(0:nx_max,l_max) & ,dfeldo(l_max),j_f(nf_max) & ,d2s(nf_max) complex*16 i,selm(l_max,-l_max:l_max),z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /excf/ selm,dfeldo common /redtrans/ dbde,dbdel common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /d2s_cd/ d2s C tmp1 = dfloat(z)/(hbarc*v) tmp1 = 10.d00*tmp1*tmp1*e2 C factor 10: fm^2 -> mb tmp2(0) = 0.d00 tmp2(1) = 1.d00 do 100 i1=2,l_max,1 tmp2(i1) = tmp2(i1-1)/(a*a) 100 continue C C +++++++++++++++++++++++++++++++++++++++ call get_dfeldo(e,v,a,theta,case,lamax) C +++++++++++++++++++++++++++++++++++++++ C do 110 i1=1,nf,1 if (incl(i1).eq.1) then i2 = lambda_f(i1) if (el(i2).eq.1) then d2s(i1) = tmp1*tmp2(i2)*dbde(ne,i1)*dfeldo(i2) endif else d2s(i1) = 0.d00 endif 110 continue C return end C ******************************************** subroutine get_sigma2(nf,nx,e_i,j_i,j_b,j_c) C ******************************************** implicit none integer*4 nf,nx,nx_max,nf_max,l_max,n_max parameter (nx_max=501,nf_max=20,l_max=3,n_max=5000) integer*4 l_f(nf_max),incl(nf_max),el(0:l_max) & ,i1,i2,i3,i4,lambda_f(nf_max) real*8 mu_bc,e_i,j_i,j_b,j_c,j_f(nf_max) & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu,arg & ,dbde(0:nx_max,nf_max),dbdel(0:nx_max,l_max) & ,h2m,zze2,k_g,k_bc,tmp1,tmp2,tmp3,tmp4,fac(l_max) & ,xv(0:nx_max),wf_f(0:n_max) & ,sig_ph(0:nx_max,nf_max) & ,sig_cp(0:nx_max,nf_max) & ,sig_phl(0:nx_max,0:l_max),sig_cpl(0:nx_max,0:l_max) & ,sig_pht(0:nx_max),sig_cpt(0:nx_max) & ,sfac(0:nx_max),sfacl(0:nx_max,0:l_max) complex*16 i,wffac,z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /wave_f/ xv,wf_f,wffac common /redtrans/ dbde,dbdel common /sigma2/ sig_ph,sig_cp,sig_phl,sig_cpl & ,sig_pht,sig_cpt,sfac,sfacl common /par_bc/ mu_bc,h2m,zze2 C C factor 10: fm^2 -> mb tmp1 = tpi*tpi*tpi*10.d00 C do 100 i1=1,l_max,1 if (el(i1).eq.1) then fac(i1) = 1.d00 do 110 i2=(2*i1+1),1,(-2) fac(i1) = fac(i1)*dfloat(i2*i2) 110 continue fac(i1) = tmp1*dfloat(i1+1)/(fac(i1)*dfloat(i1)) endif c write(6,*) i1,fac(i1) 100 continue C do 200 i1=1,nf,1 if (incl(i1).eq.1) then i3 = lambda_f(i1) C??? ccc tmp3 = 2.d00*(2.d00*j_i+1.d00)/(2.d00*j_f(i1)+1.d00) C??? tmp3 = 2.d00*(2.d00*j_i+1.d00) & /((2.d00*j_b+1.d00)*(2.d00*j_c+1.d00)) C do 210 i2=0,nx,1 if (xv(i2).le.0.d00) then sig_ph(i2,i1) = 0.d00 sig_cp(i2,i1) = 0.d00 else k_g = (e_i+xv(i2))/hbarc k_bc = dsqrt(xv(i2)/h2m) tmp2 = fac(i3) do 220 i4=1,(2*i3-1),1 tmp2 =tmp2*k_g 220 continue sig_ph(i2,i1) = tmp2*dbde(i2,i1) tmp4 = tmp3*k_g*k_g/(k_bc*k_bc) sig_cp(i2,i1) = sig_ph(i2,i1)*tmp4 endif 210 continue endif 200 continue C do 300 i1=0,nx,1 do 320 i2=1,l_max,1 sig_phl(i1,i2) = 0.d00 sig_cpl(i1,i2) = 0.d00 320 continue do 310 i2=1,nf,1 if (incl(i2).eq.1) then sig_phl(i1,lambda_f(i2)) = sig_phl(i1,lambda_f(i2)) & +sig_ph(i1,i2) sig_cpl(i1,lambda_f(i2)) = sig_cpl(i1,lambda_f(i2)) & +sig_cp(i1,i2) endif 310 continue 300 continue C do 400 i1=0,nx,1 sig_pht(i1) = 0.d00 sig_cpt(i1) = 0.d00 do 410 i2=1,l_max sig_pht(i1) = sig_pht(i1)+sig_phl(i1,i2) sig_cpt(i1) = sig_cpt(i1)+sig_cpl(i1,i2) 410 continue 400 continue C do 500 i1=0,nx,1 if (xv(i1).gt.0.d00) then arg = hbarc*dsqrt(2.d00*xv(i1)/mu_bc) arg = tpi*zze2/arg arg = xv(i1)*dexp(arg) sfac(i1) = sig_cpt(i1)*arg do 510 i2=1,l_max,1 sfacl(i1,i2) = sig_cpl(i1,i2)*arg 510 continue else sfac(i1) = 0.d00 do 520 i2=1,l_max,1 sfacl(i1,i2) = sig_cpl(i1,i2)*arg 520 continue endif 500 continue C if (xv(0).eq.0.d00) then if (nx.ge.4) then sfac(0) = 4.d00*(sfac(1)+sfac(3))-6.d00*sfac(2)-sfac(4) if (sfac(0).lt.0.d00) sfac(0) = 0.d00 do 530 i1=1,l_max,1 sfacl(0,i1) = 4.d00*(sfacl(1,i1)+sfacl(3,i1)) & -6.d00*sfacl(2,i1)-sfacl(4,i1) if (sfacl(0,i1).lt.0.d00) sfacl(0,i1) = 0.d00 530 continue endif endif C return end C ************************************ subroutine get_dbde(nf,nx,mu_bc,j_i) C ************************************ implicit none integer*4 nf,nx,nx_max,nf_max,l_max,n_max parameter (nx_max=501,nf_max=20,l_max=3,n_max=5000) integer*4 l_f(nf_max),incl(nf_max),el(0:l_max) & ,i1,i2,i3,lambda_f(nf_max) real*8 mu_bc,j_i,j_f(nf_max) & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,xv(0:nx_max),wf_f(0:n_max),tmp1,tmp2,tmp3 & ,dbde(0:nx_max,nf_max),p_bc & ,dbdel(0:nx_max,l_max) complex*16 i,rmat(0:nx_max,nf_max),wffac,z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /wave_f/ xv,wf_f,wffac common /redmat/ rmat common /redtrans/ dbde,dbdel C tmp1 = tpi*hbarc tmp1 = mu_bc/(tmp1*tmp1*tmp1) C do 200 i1=0,nx,1 do 210 i2=1,nf,1 dbde(i1,i2) = 0.d00 210 continue 200 continue C do 100 i1=0,nx,1 p_bc = dsqrt(2.d00*mu_bc*xv(i1)) tmp2 = tmp1*p_bc do 110 i2=1,nf,1 if (incl(i2).eq.1) then tmp3 = tmp2*(2.d00*j_f(i2)+1.d00)/(2.d00*j_i+1.d00) i3 = lambda_f(i2) if (el(i3).eq.1) then dbde(i1,i2) = tmp3 & *dreal(dconjg(rmat(i1,i2))*rmat(i1,i2)) else dbde(i1,i2) = 0.d00 endif endif 110 continue 100 continue C do 320 i3=1,l_max,1 do 330 i1=0,nx,1 dbdel(i1,i3) = 0.d00 330 continue 320 continue do 300 i2=1,nf,1 if (incl(i2).eq.1) then i3 = lambda_f(i2) if (el(i3).eq.1) then do 310 i1=0,nx,1 dbdel(i1,i3) = dbdel(i1,i3)+dbde(i1,i2) 310 continue endif endif 300 continue C return end C ********************************** subroutine get_wf_scatt(ne,c,nr,s) C ********************************** implicit none integer*4 n_max,nx_max,l_max,nf_max parameter (n_max=5000,nx_max=501,l_max=3,nf_max=20) integer*4 ne,c,nr,l_f(nf_max),incl(nf_max),el(0:l_max),nm & ,z_b,z_c,i1,ll,iexp(30),lambda_f(nf_max) real*8 er,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,hi1,hi2,h2,h,twh,twh2,rv(l_max,0:n_max) & ,xv(0:nx_max),wf_f(0:n_max),j_f(nf_max) & ,zze2,k,s,dll,ls,rc,ac,vc,rso,aso,fso,h2m,mu_bc & ,vc_f(nf_max),rrc_f(nf_max),ac_f(nf_max) & ,fso_f(nf_max),rrso_f(nf_max),aso_f(nf_max) & ,tmp1(0:n_max),tmp2(0:n_max),pot(0:n_max),tmp3 & ,logder,kmat,delta(nf_max,0:nx_max) & ,ddelta(nf_max,0:nx_max) & ,v,eta,arg,f(30),fp(30),g(30),gp(30),sigma(30) complex*16 i,smat,wffac,z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /wave_f/ xv,wf_f,wffac common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /pot_f/ vc_f,rrc_f,ac_f,fso_f,rrso_f,aso_f common /effpot/ pot common /par_bc/ mu_bc,h2m,zze2 common /phases/ delta,ddelta C er = xv(ne) if (er.le.0.d00) then delta(c,ne) = 0.d00 wffac = 0.d00 do 200 i1=0,nr,1 wf_f(i1) = 0.d00 200 continue else k = dsqrt(er/h2m) dll = dfloat(l_f(c)*(l_f(c)+1)) ls = (j_f(c)*(j_f(c)+1.d00)-dfloat(l_f(c)*(l_f(c)+1)) & -s*(s+1.d00))/2.d00 C rc = rrc_f(c) ac = ac_f(c) vc = vc_f(c) rso = rrso_f(c) aso = aso_f(c) fso = fso_f(c) C nm = idnint((rc+10.d00)/h) if (nm.gt.nr) nm = nr-3 C C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ C do 100 i1=1,(nr-1),1 tmp1(i1) = pot(i1)+er/h2m 100 continue C tmp2(0) = 0.d00 tmp2(1) = 1.d00 do 110 i1=1,(nm+1),1 tmp2(i1+1) = (tmp2(i1)*(hi2-10.d00*tmp1(i1)) & -tmp2(i1-1)*(hi1+tmp1(i1-1))) & /(hi1+tmp1(i1+1)) 110 continue C logder = (8.d00*(tmp2(nm+1)-tmp2(nm-1)) & -tmp2(nm+2)+tmp2(nm-2))/(twh*tmp2(nm)) C v = dsqrt(2.d00*er/mu_bc) eta = zze2/(hbarc*v) ll = l_f(c) arg = k*h*dfloat(nm) C ++++++++++++++++++++++++++++++++++++++++++++ call dfcoul(ll,eta,arg,f,fp,g,gp,iexp,sigma) C ++++++++++++++++++++++++++++++++++++++++++++ C kmat = (logder*f(ll+1)-k*fp(ll+1))/(k*gp(ll+1)-logder*g(ll+1)) delta(c,ne) = datan(kmat) smat = cdexp(2.d00*i*delta(c,ne)) delta(c,ne) = delta(c,ne)/radfac if (delta(c,ne).lt.0.d00) delta(c,ne) = delta(c,ne)+180.d00 C wffac = fpi*cdexp(i*sigma(ll+1))*(smat+1.d00)/(2.d00*k) do 120 i1=1,ll,1 wffac= wffac/i 120 continue C tmp3 = (f(ll+1)+kmat*g(ll+1))/tmp2(nm) do 130 i1=0,nm,1 wf_f(i1) = tmp3*tmp2(i1) 130 continue do 140 i1=(nm+1),nr,1 arg = k*dfloat(i1)*h C ++++++++++++++++++++++++++++++++++++++++++++ call dfcoul(ll,eta,arg,f,fp,g,gp,iexp,sigma) C ++++++++++++++++++++++++++++++++++++++++++++ wf_f(i1) = f(ll+1)+kmat*g(ll+1) 140 continue endif C return end C **************************************** subroutine get_redmat(n_i,nf,nx,nr,s,de) C **************************************** implicit none integer*4 l_max,n_max,nb_max,nf_max,nx_max parameter (l_max=3,n_max=5000,nb_max=5,nf_max=20,nx_max=501) integer*4 i1,i2,i3,i4,n_i,nr,nf,nx,z_b,z_c & ,l_f(nf_max),incl(nf_max),el(0:l_max),lambda_f(nf_max) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu,h2m,zze2 & ,res,tmp(0:n_max),er,mu_bc,s,j_f(nf_max),de & ,hi1,hi2,h2,h,twh,twh2,rv(l_max,0:n_max) & ,wf_i(nb_max,0:n_max),eb(nb_max) & ,xv(0:nx_max),wf_f(0:n_max),coeff_d(l_max,nf_max) complex*16 i,rmat(0:nx_max,nf_max),wffac,z_eff(l_max),tmp1 C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /wave_i/ eb,wf_i common /wave_f/ xv,wf_f,wffac common /co_d/ coeff_d common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /redmat/ rmat common /par_bc/ mu_bc,h2m,zze2 C do 200 i1=0,nx,1 do 210 i2=1,nf,1 rmat(i1,i2) = 0.d00 210 continue 200 continue C do 100 i1=0,nx,1 if (xv(i1).gt.0.d00) then do 110 i2=1,nf,1 if (incl(i2).eq.1) then C +++++++++++++++++++++++++++++ call get_wf_scatt(i1,i2,nr,s) C +++++++++++++++++++++++++++++ c write(6,*) 'wf_scatt calculated' i3 = lambda_f(i2) if (el(i3).eq.1) then tmp1 = dsqrt(e2)*z_eff(i3) c write(6,*) 'er',xv(i1) do 130 i4=0,nr,1 tmp(i4) = rv(i3,i4)*wf_i(n_i,i4)*wf_f(i4) c write(31,*) rv(1,i4),wf_f(i4),tmp(i4) 130 continue c stop C ++++++++++++++++++++++++++ call simpson(nr,tmp,h,res) C ++++++++++++++++++++++++++ rmat(i1,i2) = res*tmp1*wffac*coeff_d(i3,i2) c rmat(i1,i2) = 1.d00 endif endif 110 continue endif 100 continue C return end C ************************************ subroutine get_coeff_b(nf,l_i,s,j_i) C ************************************ implicit none integer*4 nf,l_i,l_max,nf_max parameter (l_max=3,nf_max=20) integer*4 l_f(nf_max),incl(nf_max),el(0:l_max) & ,i1,i2,i3,i4,i5,lambda_f(nf_max) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,s,j_i,j_f(nf_max),dl_i,dl_f,dl_fp,dj_f,dj_fp & ,dla,dlap,dl,clego,sjsym,tmp & ,coeff_b(nf_max,nf_max,0:(2*l_max)) complex*16 i,z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /co_b/ coeff_b C do 100 i1=1,nf,1 do 110 i2=1,nf,1 do 120 i5=0,(2*l_max),1 coeff_b(i1,i2,i5) = 0.d00 120 continue 110 continue 100 continue C dl_i = dfloat(l_i) C do 200 i1=1,nf,1 if (incl(i1).eq.1) then dl_f = dfloat(l_f(i1)) dj_f = j_f(i1) i3 = lambda_f(i1) if (el(i3).eq.1) then dla = dfloat(i3) do 210 i2=1,nf,1 if (incl(i2).eq.1) then dl_fp = dfloat(l_f(i2)) dj_fp = j_f(i2) i4 = lambda_f(i2) if (el(i4).eq.1) then dlap = dfloat(i4) do 240 i5=0,(2*l_max),1 dl = dfloat(i5) tmp = (2.d00*dl_f+1.d00)*(2.d00*dl_fp+1.d00) & /(fpi*(2.d00*dla+1.d00)) tmp = dsqrt(tmp)*(2.d00*dj_f+1.d00)*(2.d00*dj_fp+1.d00) & /((2.d00*dla+1.d00)*(2.d00*dlap+1.d00)) tmp = tmp*clego(dl_f,0.d00,dl_fp,0.d00,dl,0.d00) & *sjsym(dl_fp,s,dj_fp,dj_f,dl,dl_f) & *sjsym(dj_fp,j_i,dlap,dla,dl,dj_f) if (mod(idnint(s+dj_fp+dj_f+j_i+dlap),2).ne.0) then coeff_b(i1,i2,i5) = -tmp/(2.d00*j_i+1.d00) else coeff_b(i1,i2,i5) = tmp/(2.d00*j_i+1.d00) endif 240 continue endif endif 210 continue endif endif 200 continue C return end C ************************************ subroutine get_coeff_d(nf,l_i,s,j_i) C ************************************ implicit none integer*4 nf,l_i,l_max,nf_max parameter (l_max=3,nf_max=20) integer*4 el(0:l_max),incl(nf_max),l_f(nf_max) & ,lambda_f(nf_max),i1,i2 real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,s,j_i,j_f(nf_max),dl_i,dl,dl_f,dj_f & ,tmp1,tmp2(l_max),coeff_d(l_max,nf_max) & ,clego,sjsym complex*16 i,z_eff(l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /fstates/ z_eff,j_f,l_f,lambda_f,incl,el common /co_d/ coeff_d C dl_i = dfloat(l_i) tmp1 = dsqrt((2.d00*dl_i+1.d00)*(2.d00*j_i+1.d00)) if (mod(idnint(s+dl_i+j_i),2).ne.0) tmp1 = -tmp1 C do 100 i1=1,l_max,1 if (el(i1).eq.1) then tmp2(i1) = dsqrt(dfloat(2*i1+1)/fpi) else tmp2(i1) = 0.d00 endif 100 continue C do 110 i1=1,l_max,1 dl = dfloat(i1) do 120 i2=1,nf,1 dl_f = dfloat(l_f(i2)) dj_f = j_f(i2) if ((el(i1).eq.1).and.(incl(i2).eq.1)) then coeff_d(i1,i2) = tmp1*tmp2(i1) & *clego(dl_i,0.d00,dl,0.d00,dl_f,0.d00) & *sjsym(dl_i,s,j_i,dj_f,dl,dl_f) c write(6,*) 'i1,i2,coeff_d',i1,i2,coeff_d(i1,i2) else coeff_d(i1,i2) = 0.d00 endif 120 continue 110 continue C return end C ******************************************************* subroutine get_wf_bound(nr,z_b,z_c,nb,l,s,j,prec & ,rc,ac,vc,rso,aso,fso,e,case,nc) C ******************************************************* implicit none integer*4 nr,z_b,z_c,n,l,case,l_max,n_max,nb_max,nc_max parameter (l_max=3,n_max=5000,nb_max=5,nc_max=2000) integer*4 i1,i2,nm,nc,nb,ns,nk real*8 s,j,mu_bc,rc,ac,vc,rso,aso,fso,e & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,zze2,h2m,dll,ls,hi1,hi2,h2,h,twh,twh2 & ,eb_old,de,prec,dld & ,rv(l_max,0:n_max),wf_i(nb_max,0:n_max) & ,eb(nb_max),vcv(0:2),fv(0:2),df & ,pot(0:n_max) complex*16 i C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /wave_i/ eb,wf_i common /par_bc/ mu_bc,h2m,zze2 C common /effpot/ pot C dll = dfloat(l*(l+1)) ls = (j*(j+1.d00)-dfloat(l*(l+1))-s*(s+1.d00))/2.d00 nm = idnint(rc/h) C if (nb.ne.0) then C ++++++++++++++++++++ call init_wf(nr,rc) C ++++++++++++++++++++ if (case.eq.1) then C ***** calculate binding energy from given potential ***** C by variational + shooting method C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ do 140 i2=1,nb,1 nc = 0 de = 1.d99 120 continue nc = nc+1 eb_old = eb(i2) do 100 i1=1,50,1 C +++++++++++++++++++++ call get_wf_sw(i2,nr) C +++++++++++++++++++++ 100 continue de = dabs(eb(i2)-eb_old) if (nc.gt.nc_max) goto 130 if (de.gt.(1.d03*prec)) goto 120 nc = 0 130 continue 140 continue C vcv(0) = eb(nb) C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ fv(0) = dld eb(nb) = eb(nb)-0.1d00 vcv(1) = eb(nb) C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ fv(1) = dld 110 continue vcv(2) = (vcv(1)*fv(0)-vcv(0)*fv(1))/(fv(0)-fv(1)) eb(nb) = vcv(2) C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ fv(2) = dld C vcv(0) = vcv(1) vcv(1) = vcv(2) fv(0) = fv(1) fv(1) = fv(2) if (dabs(dld).gt.1.d-12) goto 110 eb(nb) = vcv(1) C e = -eb(nb) else C ***** calculate potential depth for given energy ***** C by shooting method eb(nb) = -e vcv(0) = 50.d00 210 continue vc = vcv(0) C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ c write(6,*) vcv(0),nk c stop fv(0) = dld if (nk.lt.nb) then vcv(0) = vcv(0)+0.1d00 goto 210 endif if (nk.gt.nb) then vcv(0) = vcv(0)-0.1d00 goto 210 endif C vcv(1) = vcv(0)+0.1d00 vc = vcv(1) C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ fv(1) = dld C 200 continue vcv(2) = (vcv(1)*fv(0)-vcv(0)*fv(1))/(fv(0)-fv(1)) vc = vcv(2) C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ fv(2) = dld C vcv(0) = vcv(1) vcv(1) = vcv(2) fv(0) = fv(1) fv(1) = fv(2) df = dabs(fv(1)-fv(0)) if (df.gt.1.d-12) goto 200 C endif else nc = -1 endif C c do 600 i1=0,nr,1 c write(42,*) rv(1,i1),wf_i(nb,i1) c write(33,*) rv(1,i1),pot(i1) c 600 continue c stop C return end C ************************************* subroutine get_wf_sh(nm,nb,nr,dld,nk) C ************************************* implicit none integer*4 nr,i1,i2,n_max,l_max,nb_max,nb,nm,nk parameter (n_max=5000,l_max=3,nb_max=5) real*8 h2m,hi1,hi2,h2,h,twh,twh2,res,fac,ld2,ld3,dld & ,rv(l_max,0:n_max),pot(0:n_max),er,mu_bc,zze2 & ,eb(nb_max),wf_i(nb_max,0:n_max) & ,tmp1(0:n_max),tmp2(0:n_max),tmp3(0:n_max) C common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /effpot/ pot common /wave_i/ eb,wf_i common /par_bc/ mu_bc,h2m,zze2 C ccc if (dabs(wf_i(nb,nr-1)).lt.1.d-30) return C C 10 keV Grenze fuer Bindingsenergie ??? c if (eb(nb).gt.(-0.01d00)) then er = eb(nb) c else c er = -0.01d00 c endif C c write(6,*) 'er',er do 100 i1=1,(nr-1),1 tmp1(i1) = pot(i1)+er/h2m 100 continue C tmp2(0) = 0.d00 tmp2(1) = 1.d00 cccwf_i(nb,1) do 110 i1=1,(nm+2),1 tmp2(i1+1) = (tmp2(i1)*(hi2-10.d00*tmp1(i1)) & -tmp2(i1-1)*(hi1+tmp1(i1-1))) & /(hi1+tmp1(i1+1)) 110 continue tmp3(nr) = 0.d00 tmp3(nr-1) = 1.d00 cccwf_i(nb,nr-1) do 120 i1=(nr-1),(nm-2),(-1) tmp3(i1-1) = (tmp3(i1)*(hi2-10.d00*tmp1(i1)) & -tmp3(i1+1)*(hi1+tmp1(i1+1))) & /(hi1+tmp1(i1-1)) 120 continue C if (tmp3(nm).ne.0.d00) then fac = tmp2(nm)/tmp3(nm) else fac = tmp2(nm+1)/tmp3(nm+1) endif C if (tmp2(nm)*tmp3(nm).ne.0.d00) then ld2 = (8.d00*(tmp2(nm+1)-tmp2(nm-1)) & -tmp2(nm+2)+tmp2(nm-2))/(twh*tmp2(nm)) ld3 = (8.d00*(tmp3(nm+1)-tmp3(nm-1)) & -tmp3(nm+2)+tmp3(nm-2))/(twh*tmp3(nm)) else ld2 = (8.d00*(tmp2(nm+2)-tmp2(nm)) & -tmp2(nm+3)+tmp2(nm-1))/(twh*tmp2(nm+1)) ld3 = (8.d00*(tmp3(nm+2)-tmp3(nm)) & -tmp3(nm+3)+tmp3(nm-1))/(twh*tmp3(nm+1)) endif dld = ld2-ld3 C do 130 i1=0,(nm-1),1 wf_i(nb,i1) = tmp2(i1) tmp2(i1) = tmp2(i1)*tmp2(i1) 130 continue do 140 i1=nm,nr,1 wf_i(nb,i1) = fac*tmp3(i1) tmp2(i1) = wf_i(nb,i1)*wf_i(nb,i1) 140 continue C C +++++++++++++++++++++++++++ call simpson(nr,tmp2,h,res) C +++++++++++++++++++++++++++ c write(6,*) 'res',res res = 1.d00/dsqrt(res) nk = 1 do 150 i1=1,(nr-1),1 wf_i(nb,i1) = wf_i(nb,i1)*res if (wf_i(nb,i1)*wf_i(nb,i1+1).lt.0.d00) nk = nk+1 c write(41,*) rv(1,i1),wf_i(nb,i1) 150 continue C return end C *************************** subroutine get_wf_sw(nb,nr) C *************************** implicit none integer*4 nr,i1,i2,n_max,l_max,nb_max,nb parameter (n_max=5000,l_max=3,nb_max=5) real*8 h2m,hi1,hi2,h2,h,twh,twh2,res,relax_wf,relax_eb & ,rv(l_max,0:n_max),pot(0:n_max),mu_bc,zze2 & ,eb(nb_max),wf_i(nb_max,0:n_max) & ,tmp1(0:n_max),tmp2(0:n_max) C common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /effpot/ pot common /wave_i/ eb,wf_i common /par_bc/ mu_bc,h2m,zze2 C relax_wf = 0.5d00 c relax_eb = 1.0d00 C do 100 i1=1,(nr-1),1 tmp1(i1) = pot(i1)+eb(nb)/h2m 100 continue c stop C do 110 i2=1,3,1 do 120 i1=1,(nr-1),1 tmp2(i1) = (wf_i(nb,i1+1)*(hi1+tmp1(i1-1)) & +wf_i(nb,i1-1)*(hi1+tmp1(i1+1))) & /(hi2-10.d00*tmp1(i1)) 120 continue do 130 i1=1,(nr-1),1 wf_i(nb,i1) = relax_wf*tmp2(i1) & +(1.d00-relax_wf)*wf_i(nb,i1) 130 continue 110 continue C do 200 i2=1,(nb-1),1 do 210 i1=0,nr,1 tmp2(i1) = wf_i(nb,i1)*wf_i(i2,i1) 210 continue C +++++++++++++++++++++++++++ call simpson(nr,tmp2,h,res) C +++++++++++++++++++++++++++ do 220 i1=0,nr,1 wf_i(nb,i1) = wf_i(nb,i1)-res*wf_i(i2,i1) 220 continue 200 continue C do 140 i1=0,nr,1 tmp2(i1) = wf_i(nb,i1)*wf_i(nb,i1) 140 continue C +++++++++++++++++++++++++++ call simpson(nr,tmp2,h,res) C +++++++++++++++++++++++++++ res = 1.d00/dsqrt(res) do 150 i1=0,nr,1 wf_i(nb,i1) = wf_i(nb,i1)*res 150 continue C i1 = 1 tmp2(i1) = (11.d00*wf_i(nb,i1-1)-20.d00*wf_i(nb,i1) & +6.d00*wf_i(nb,i1+1)+4.d00*wf_i(nb,i1+2) & -wf_i(nb,i1+3))/twh2 tmp2(i1) = (tmp2(i1)+pot(i1)*wf_i(nb,i1)) & *wf_i(nb,i1) do 160 i1=2,(nr-2),1 c tmp2(i1) = (wf_i(nb,i1+1)+wf_i(nb,i1-1) c & -2.d00*wf_i(nb,i1))/h2 tmp2(i1) = (16.d00*(wf_i(nb,i1+1)+wf_i(nb,i1-1)) & -wf_i(nb,i1+2)-wf_i(nb,i1+2) & -30.d00*wf_i(nb,i1))/twh2 tmp2(i1) = (tmp2(i1)+pot(i1)*wf_i(nb,i1)) & *wf_i(nb,i1) 160 continue i1 = nr-1 tmp2(i1) = (11.d00*wf_i(nb,i1+1)-20.d00*wf_i(nb,i1) & +6.d00*wf_i(nb,i1-1)+4.d00*wf_i(nb,i1-2) & -wf_i(nb,i1-3))/twh2 tmp2(i1) = (tmp2(i1)+pot(i1)*wf_i(nb,i1)) & *wf_i(nb,i1) C +++++++++++++++++++++++++++ call simpson(nr,tmp2,h,res) C +++++++++++++++++++++++++++ eb(nb) = -res*h2m ccc eb(nb) = -relax_eb*res*h2m+(1.d00-relax_eb)*eb(nb) c write(6,*) eb(nb) C return end C ************************* subroutine init_wf(nr,rc) C ************************* implicit none integer*4 nr,n_max,l_max,nb_max,i1,i2 parameter (n_max=5000,l_max=3,nb_max=5) real*8 rc,hi1,hi2,h2,h,twh,twh2,fac,k,res & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,rv(l_max,0:n_max),eb(nb_max) & ,tmp1(0:n_max),wf_i(nb_max,0:n_max) & ,tmp2(0:n_max) complex*16 i C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /wave_i/ eb,wf_i C fac = pi/(nr*h) do 130 i1=1,nr,1 tmp2(i1) = dexp(-rv(1,i1)/rc) ccc tmp2(i1) = dexp(-rv(1,i1)/(10.d00*rv(1,nr))) 130 continue do 100 i1=1,nb_max,1 eb(i1) = 0.d00 k = dfloat(i1)*fac do 110 i2=0,nr,1 wf_i(i1,i2) = dsin(k*rv(1,i2))*tmp2(i2) tmp1(i2) = wf_i(i1,i2)*wf_i(i1,i2) 110 continue C +++++++++++++++++++++++++++ call simpson(nr,tmp1,h,res) C +++++++++++++++++++++++++++ res = 1.d00/dsqrt(res) do 120 i2=0,nr,1 wf_i(i1,i2) = res*wf_i(i1,i2) 120 continue 100 continue C return end C ***************************************************** subroutine get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C ***************************************************** implicit none integer*4 n_max,l_max,nr,i1 parameter (n_max=5000,l_max=3) real*8 dll,zze2,h2m,ls,rc,ac,vc,rso,aso,fso,mu_bc & ,rv(l_max,0:n_max),pot(0:n_max) & ,fws,r,rr,a,fac,tmp,x,hi1,hi2,h2,h,twh,twh2 C common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /effpot/ pot common /par_bc/ mu_bc,h2m,zze2 C fws(r,rr,a) = 1.d00/(1.d00+dexp((r-rr)/a)) C fac = fso*ls/aso do 100 i1=1,nr,1 r = rv(1,i1) tmp = fws(r,rso,aso) tmp = vc*(fws(r,rc,ac)+fac*tmp*tmp*dexp((r-rso)/aso)/r) if (r.gt.rc) then tmp = tmp-zze2/r else x = r/rc tmp = tmp-zze2*(3.d00-x*x)/(2.d00*rc) endif pot(i1) = tmp/h2m-dll/rv(2,i1) ccc pot(i1) = 0.d00 c write(30,*) r,(-pot(i1)*h2m) 100 continue C return end C ********************************** subroutine kinematics(e,a,v,gamma) C ********************************** implicit none integer*4 z_x,z_a real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,m_x,m_a,mu_ax & ,t_p,e,vi,vf,v,a,gamma & ,s,e_ax,p_ax,m_aa,mu_aax,e_aax,p_aax complex *16 i C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /par_ax/ z_x,z_a,m_x,m_a,mu_ax,t_p C s = (m_x+m_a)*(m_x+m_a)+2.d00*m_x*t_p e_ax = dsqrt(s)-m_a-m_x+mu_ax p_ax = dsqrt(e_ax*e_ax-mu_ax*mu_ax) vi = p_ax/e_ax m_aa = m_a+e mu_aax = m_aa*m_x/(m_aa+m_x) e_aax = dsqrt(s)-m_aa-m_x+mu_aax p_aax = dsqrt(e_aax*e_aax-mu_aax*mu_aax) vf = p_aax/e_aax C v = dsqrt(vi*vf) gamma = 1.d00/dsqrt(1.d00-v*v) a = dfloat(z_x*z_a)*e2/(mu_ax*v*v) C return end C ****************************** subroutine err_msg(file_out,n) C ****************************** implicit none integer*4 n,file_out C write(file_out,*) '************ error ************' if (n.eq.1) write(file_out,*) 'nf > nf_max' if (n.eq.2) write(file_out,*) 'nth > nth_max' if (n.eq.3) write(file_out,*) 'ne > nx_max' if (n.eq.4) write(file_out,*) '2*np+1 > nx_max' if (n.eq.5) write(file_out,*) 'nth_max > n_max' if (n.eq.6) write(file_out,*) 'no bound state found' c if (n.eq.7) write(file_out,*) '2*np_max+1 > nx_max' if (n.eq.8) write(file_out,*) 'nr > n_max' if (n.eq.9) write(file_out,*) & 'bound state wave function not precise' if (n.eq.10) write(file_out,*) & 'principal quantum number too small' write(file_out,*) '***** calculation stopped *****' C stop end C ********************* subroutine init_const C ********************* implicit none real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu complex*16 i C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu C i = (0.d00,1.d00) C pi = 4.d00*datan(1.d00) tpi = 2.d00*pi fpi = 4.d00*pi pih = 0.5d00*pi radfac = pi/180.d00 e2 = 1.44d00 hbarc = 197.3286d00 amu = 931.49386d00 C return end C ******************************************************* subroutine get_wf_boundold(nr,z_b,z_c,nb,l,s,j,prec & ,rc,ac,vc,rso,aso,fso,e,case,nc) C ******************************************************* implicit none integer*4 nr,z_b,z_c,n,l,case,l_max,n_max,nb_max,nc_max parameter (l_max=3,n_max=5000,nb_max=5,nc_max=2000) integer*4 i1,i2,nm,nc,nb,ns,nk real*8 s,j,mu_bc,rc,ac,vc,rso,aso,fso,e & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,zze2,h2m,dll,ls,hi1,hi2,h2,h,twh,twh2 & ,eb_old,de,prec,dld & ,rv(l_max,0:n_max),wf_i(nb_max,0:n_max) & ,eb(nb_max),vcv(0:2),fv(0:2),df complex*16 i C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /mesh/ rv,hi1,hi2,h2,h,twh,twh2 common /wave_i/ eb,wf_i common /par_bc/ mu_bc,h2m,zze2 C dll = dfloat(l*(l+1)) ls = (j*(j+1.d00)-dfloat(l*(l+1))-s*(s+1.d00))/2.d00 nm = idnint(rc/h) C if (nb.ne.0) then C ++++++++++++++++++++ call init_wf(nr,rc) C ++++++++++++++++++++ if (case.eq.1) then C ***** calculate binding energy from given potential ***** C by variational + shooting method C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ do 140 i2=1,nb,1 nc = 0 ns = 0 de = 1.d99 120 continue nc = nc+1 c write(6,*) eb(i2) eb_old = eb(i2) do 100 i1=1,50,1 C +++++++++++++++++++++ call get_wf_sw(i2,nr) C +++++++++++++++++++++ 100 continue if (ns.eq.0) then if (de.lt.1.d02*prec) then ns = 1 C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,i2,nr,dld,nk) C +++++++++++++++++++++++++++++++ endif endif do 110 i1=1,50,1 C +++++++++++++++++++++ call get_wf_sw(i2,nr) C +++++++++++++++++++++ 110 continue de = dabs(eb(i2)-eb_old) if (nc.gt.nc_max) goto 130 if (de.gt.prec) goto 120 nc = 0 130 continue 140 continue e = -eb(nb) else C ***** calculate potential depth for given energy ***** C by shooting method eb(nb) = -e vcv(0) = 50.d00 210 continue vc = vcv(0) C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ c write(6,*) vcv(0),nk c stop fv(0) = dld if (nk.lt.nb) then vcv(0) = vcv(0)+1.d00 goto 210 endif if (nk.gt.nb) then vcv(0) = vcv(0)-1.d00 goto 210 endif C vcv(1) = vcv(0)-1.d00 vc = vcv(1) C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ fv(1) = dld C 200 continue vcv(2) = (vcv(1)*fv(0)-vcv(0)*fv(1))/(fv(0)-fv(1)) vc = vcv(2) C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ C +++++++++++++++++++++++++++++++ call get_wf_sh(nm,nb,nr,dld,nk) C +++++++++++++++++++++++++++++++ fv(2) = dld C vcv(0) = vcv(1) vcv(1) = vcv(2) fv(0) = fv(1) fv(1) = fv(2) df = dabs(fv(1)-fv(0)) if (df.gt.1.d-08) goto 200 C vc = vcv(1) vcv(0) = vc C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ do 500 i2=1,nb,1 nc = 0 eb(i2) = 1.d99 410 continue nc = nc+1 eb_old = eb(i2) do 320 i1=1,50,1 C +++++++++++++++++++++ call get_wf_sw(i2,nr) C +++++++++++++++++++++ 320 continue de = dabs(eb(i2)-eb_old) if (nc.gt.nc_max) goto 530 if (de.gt.100.d00*prec) goto 410 nc = 0 530 continue 500 continue fv(0) = eb(nb) C vcv(1) = vcv(0)-0.1d00 vc = vcv(1) C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ do 510 i2=1,nb,1 nc = 0 eb(i2) = 1.d99 420 continue nc = nc+1 eb_old = eb(i2) do 330 i1=1,50,1 C +++++++++++++++++++++ call get_wf_sw(i2,nr) C +++++++++++++++++++++ 330 continue de = dabs(eb(i2)-eb_old) if (nc.gt.nc_max) goto 540 if (de.gt.100.d00*prec) goto 420 nc = 0 540 continue 510 continue fv(1) = eb(nb) C 300 continue vcv(2) = (vcv(1)*(fv(0)+e) & -vcv(0)*(fv(1)+e))/(fv(0)-fv(1)) vc = vcv(2) C +++++++++++++++++++++++++++++++++++++++++++++++ call get_effpot(nr,dll,ls,rc,ac,vc,rso,aso,fso) C +++++++++++++++++++++++++++++++++++++++++++++++ do 520 i2=1,nb,1 nc = 0 eb(i2) = 1.d99 430 continue nc = nc+1 eb_old = eb(i2) do 310 i1=1,50,1 C +++++++++++++++++++++ call get_wf_sw(i2,nr) C +++++++++++++++++++++ 310 continue de = dabs(eb(i2)-eb_old) if (nc.gt.nc_max) goto 550 if (de.gt.100.d00*prec) goto 430 nc = 0 550 continue 520 continue fv(2) = eb(nb) C vcv(0) = vcv(1) vcv(1) = vcv(2) fv(0) = fv(1) fv(1) = fv(2) df = dabs(fv(1)-fv(0)) if (df.gt.1.d-03) goto 300 endif else nc = -1 endif C c do 600 i1=1,(nr-1),1 c write(32,*) rv(1,i1),wf_i(nb,i1) c 600 continue c stop C return end C ************************************** subroutine get_ylmx(theta,phi,case_tp) C ************************************** implicit none integer*4 l_max,tl_max,l,m,case_tp parameter (l_max=3,tl_max=2*l_max) real*8 theta,phi,x,fac1,fac2 & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,alp(0:tl_max,-tl_max:tl_max) complex*16 i,cfac,carg,ylm(0:tl_max,-tl_max:tl_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /ylmtp/ ylm common /alpx/ alp C x = dcos(theta*radfac) C ++++++++++++++++++++++++ call get_alpx(x,case_tp) C ++++++++++++++++++++++++ fac2 = 1.d00/dsqrt(fpi) if (case_tp.eq.0) fac2 = fac2*tpi do 10 l=0,tl_max,1 fac1 = dsqrt(dfloat(2*l+1))*fac2 ylm(l,0) = alp(l,0)*fac1 do 20 m=1,l,1 fac1 = -fac1/dsqrt(dfloat((l-m+1)*(l+m))) ylm(l,m) = fac1*alp(l,m) 20 continue 10 continue carg = i*phi*radfac fac1 = 1.d00 do 30 m=1,tl_max,1 cfac = cdexp(dfloat(m)*carg) fac1 = -fac1 do 40 l=m,tl_max,1 ylm(l, m) = ylm(l,m)*cfac ylm(l,-m) = fac1*dconjg(ylm(l,m)) 40 continue 30 continue C return end C ***************************** subroutine get_ylm(theta,phi) C ***************************** implicit none integer*4 l_max,tl_max,l,m parameter (l_max=3,tl_max=2*l_max) real*8 theta,phi,x,fac1,fac2 & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,alp(0:tl_max,-tl_max:tl_max) complex*16 i,cfac,carg,ylm(0:tl_max,-tl_max:tl_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /ylmtp/ ylm common /alpx/ alp C x = dcos(theta*radfac) C +++++++++++++++ call get_alp(x) C +++++++++++++++ fac2 = 1.d00/dsqrt(fpi) do 10 l=0,tl_max,1 fac1 = dsqrt(dfloat(2*l+1))*fac2 ylm(l,0) = alp(l,0)*fac1 do 20 m=1,l,1 fac1 = -fac1/dsqrt(dfloat((l-m+1)*(l+m))) ylm(l,m) = fac1*alp(l,m) 20 continue 10 continue carg = i*phi*radfac fac1 = 1.d00 do 30 m=1,tl_max,1 cfac = cdexp(dfloat(m)*carg) fac1 = -fac1 do 40 l=m,tl_max,1 ylm(l, m) = ylm(l,m)*cfac ylm(l,-m) = fac1*dconjg(ylm(l,m)) 40 continue 30 continue C return end C ****************************** subroutine get_alpx(x,case_tp) C ****************************** implicit none integer*4 l_max,tl_max,l,m,case_tp parameter (l_max=3,tl_max=2*l_max) real*8 x,fac,alp(0:tl_max,-tl_max:tl_max) C common /alpx/ alp C do 100 l=0,tl_max,1 do 110 m=-l,l,1 alp(l,m) = 0.d00 110 continue 100 continue C if (case_tp.eq.0) then alp(0,0) = 1.d00 alp(1,0) = x do 200 l=2,tl_max,1 alp(l,0) = (dfloat(2*l-1)*x*alp(l-1,0) & -dfloat(l-1)*alp(l-2,0))/dfloat(l) 200 continue else alp(0,0) = 2.d00 alp(1,1) = 2.d00*datan(1.d00) do 300 l=3,tl_max,2 alp(l,l) = alp(l-2,l-2)*dfloat(l*(2*l-1)*(2*l-3)) & /dfloat(l+1) 300 continue do 310 l=2,tl_max,2 alp(l,l) = alp(l-2,l-2)*dfloat(l*(2*l-1)*(2*l-3)) & /dfloat(l+1) 310 continue do 320 l=3,tl_max,1 do 330 m=l,2,(-2) alp(l,m-2) = alp(l,m)*dfloat(m-2) & /dfloat(m*(l-m+2)*(l+m-1)) 330 continue 320 continue do 340 l=1,tl_max,1 fac = 1.d00 do 350 m=1,l,1 fac = -fac/dfloat((l+m)*(l-m+1)) alp(l,-m) = alp(l,m)*fac 350 continue 340 continue endif C c do 400 l=0,tl_max,1 c do 410 m=l,(-l),(-1) c write(6,*) l,m,alp(l,m) c 410 continue c 400 continue c stop C return end C ********************* subroutine get_alp(x) C ********************* implicit none integer*4 l_max,tl_max,l,m parameter (l_max=3,tl_max=2*l_max) real*8 x,fac,alp(0:tl_max,-tl_max:tl_max) C common /alpx/ alp C alp(0,0) = 1.d00 alp(1,0) = x fac = -dsqrt(1.d00-x*x)/2.d00 do 10 l=1,tl_max,1 alp(l,-l) = alp((l-1),(-l+1))*fac/dfloat(l) 10 continue do 20 l=2,tl_max,1 alp(l,(-l+1)) = alp((l-1),(-l+2))*fac/dfloat(l-1) 20 continue C do 30 l=2,tl_max,1 do 40 m=(-l+2),0,1 alp(l,m) = (dfl