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) = (dfloat(2*l-1)*x*alp(l-1,m) & -dfloat(l-1+m)*alp(l-2,m))/dfloat(l-m) 40 continue 30 continue C do 50 l=1,tl_max,1 fac = 1.d00 do 60 m=1,l,1 fac = -fac*dfloat((l+m)*(l-m+1)) alp(l,m) = fac*alp(l,-m) 60 continue 50 continue C return end C ******************************************* subroutine get_telm(e,v,a,theta,case,lamax) C ******************************************* implicit none integer*4 case,l_max,tl_max parameter (l_max=3,tl_max=2*l_max) integer*4 i1,i2,i3,i4,i5,i6,lamax real*8 e,v,a,theta,dfeldo(l_max),clego & ,dla,dlap,dl,dmu,dmup,dm complex*16 ctmp,selm(l_max,-l_max:l_max) & ,telm(l_max,l_max,0:tl_max,-tl_max:tl_max) C common /excf/ selm,dfeldo common /excf2/ telm C do 100 i1=1,l_max,1 do 110 i2=1,l_max,1 do 120 i3=0,tl_max,1 do 130 i4=(-i3),i3,1 telm(i1,i2,i3,i4) = 0.d00 130 continue 120 continue 110 continue 100 continue C if (theta.le.0.d00) return C C +++++++++++++++++++++++++++++++++++++ call get_selm(e,v,a,theta,case,lamax) C +++++++++++++++++++++++++++++++++++++ C do 200 i1=1,l_max,1 dla = dfloat(i1) do 210 i2=1,l_max,1 dlap = dfloat(i2) do 220 i3=0,tl_max,1 if ((iabs(i1-i2).le.i3).and.(i3.le.(i1+i2))) then dl = dfloat(i3) do 230 i4=(-i3),i3,1 dm = dfloat(i4) ctmp = 0.d00 do 240 i5=(-i1),i1,1 dmu = dfloat(i5) i6 = i5-i4 if (iabs(i6).le.i2) then dmup = dfloat(i6) ctmp = ctmp+clego(dlap,dmup,dl,dm,dla,dmu) & *dconjg(selm(i2,-i6))*selm(i1,-i5) endif 240 continue if (mod((i3+i4),2).ne.0) then telm(i1,i2,i3,i4) = -ctmp else telm(i1,i2,i3,i4) = ctmp endif 230 continue endif 220 continue 210 continue 200 continue C return end C ********************************************* subroutine get_dfeldo(e,v,a,theta,case,lamax) C ********************************************* implicit none integer*4 case,l_max,i1,lamax parameter (l_max=3) real*8 e,a,v,theta,eps,xi,gamma & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,dfeldo(l_max) complex*16 i,selm(l_max,-l_max:l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /excf/ selm,dfeldo C if (theta.gt.0.d00) then if (case.eq.0) then C ***** non relativistic ***** c gamma = 1.d00/dsqrt(1.d00-v*v) xi = e*a/(hbarc*v) eps = 1.d00/dsin(theta*radfac/2.d00) C +++++++++++++++++++++++++++++++ call get_dfeldonr(xi,eps,lamax) C +++++++++++++++++++++++++++++++ else C ***** relativistic ***** C +++++++++++++++++++++++++++++++++++ call get_dfeldor(e,v,a,theta,lamax) C +++++++++++++++++++++++++++++++++++ endif else do 100 i1=1,l_max,1 dfeldo(i1) = 0.d00 100 continue endif C return end C ******************************************* subroutine get_selm(e,v,a,theta,case,lamax) C ******************************************* implicit none integer*4 case,l_max,lamax parameter (l_max=3) real*8 e,a,v,theta,eps,xi,gamma,b,k & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,dfeldo(l_max) complex*16 i,selm(l_max,-l_max:l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /excf/ selm,dfeldo C eps = 1.d00/dsin(theta*radfac/2.d00) C if (case.eq.0) then C ***** non relativistic ***** xi = e*a/(hbarc*v) C ++++++++++++++++++++ call get_drot(theta) C ++++++++++++++++++++ C +++++++++++++++++++++++++++++++++ call get_selmnr(xi,eps,v,a,lamax) C +++++++++++++++++++++++++++++++++ else C ***** relativistic ***** gamma = 1.d00/dsqrt(1.d00-v*v) b = a*(dsqrt(eps*eps-1.d00)+pih/gamma) k = e/hbarc xi = k*b/(gamma*v) C ++++++++++++++++++++++++++++++++++ call get_selmr(xi,k,v,gamma,lamax) C ++++++++++++++++++++++++++++++++++ endif C return end C **************************************** subroutine get_selmr(xi,k,v,gamma,lamax) C **************************************** implicit none integer*4 l_max,ll_max,l,m,lamax parameter (l_max=3,ll_max=l_max+1) real*8 xi,k,v,gamma,vi,fac,lfac,tfac & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,dfeldo(l_max),kn(0:l_max) & ,lpol(0:ll_max,-ll_max:ll_max) complex*16 i,selm(l_max,-l_max:l_max) & ,gelm(1:l_max,-l_max:l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /excf/ selm,dfeldo common /bessel/ kn common /rexcf/ gelm,lpol C C +++++++++++++++++++++ call get_besselkn(xi) C +++++++++++++++++++++ vi = 1.d00/v C +++++++++++++++++ call get_gelm(vi) C +++++++++++++++++ C C negative signs because of projectile excitation fac = -vi/(fpi*gamma) do 100 l=1,l_max,1 fac = -k*fac lfac = dfloat(2*l+1) tfac = fac*lfac*dsqrt(lfac) if (l.gt.lamax) then do 120 m=-l,l,1 selm(l,m) = 0.d00 120 continue else do 110 m=-l,l,1 selm(l,m) = tfac*kn(iabs(m))*gelm(l,m) 110 continue endif 100 continue C return end C ********************** subroutine get_gelm(x) C ********************** implicit none integer*4 l_max,l,m,ll_max parameter (l_max=3,ll_max=l_max+1) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,x,fac,tmp1,tmp2,lpol(0:ll_max,-ll_max:ll_max) complex*16 i,ctmp1,ctmp2,cfac,gelm(1:l_max,-l_max:l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /rexcf/ gelm,lpol C C ++++++++++++++++++ call get_legpol(x) C ++++++++++++++++++ C fac = 4.d00*dsqrt(pi)/dsqrt(x*x-1.d00) C ctmp1 = 1.d00 tmp1 = fac do 100 l=1,l_max,1 ctmp1 = i*ctmp1 tmp1 = tmp1/dfloat(2*l+1) cfac = tmp1*ctmp1/dfloat(l) ctmp2 = cfac/i tmp2 = 1.d00 do 110 m=0,(l-1),1 ctmp2 = i*ctmp2 gelm(l,m) = ctmp2*tmp2 & *(dfloat((l+1)*(l+m))*lpol(l-1,m) & -dfloat(l*(l-m+1))*lpol(l+1,m))/dfloat(2*l+1) tmp2 = tmp2/dsqrt(dfloat((l-m)*(l+m+1))) 110 continue m = l ctmp2 = i*ctmp2 c tmp2 = tmp2/dsqrt(dfloat(l+m)) c tmp2 = tmp2/dsqrt(dfloat(l+m+1)) gelm(l,m) = ctmp2*tmp2 & *(dfloat((l+1)*(l+m))*lpol(l-1,m) & -dfloat(l*(l-m+1))*lpol(l+1,m))/dfloat(2*l+1) 100 continue C do 120 l=1,l_max,1 do 130 m=1,l,2 gelm(l,-m) = -gelm(l,m) 130 continue do 140 m=2,l,2 gelm(l,-m) = gelm(l,m) 140 continue 120 continue C c do 150 l=1,l_max,1 c do 160 m=(-l),l,1 c write(6,*) l,m,gelm(l,m) c 160 continue c 150 continue C return end C ************************ subroutine get_legpol(x) C ************************ implicit none integer*4 l_max,l,m,ll_max parameter (l_max=3,ll_max=l_max+1) real*8 x,fac,lpol(0:ll_max,-ll_max:ll_max) complex*16 gelm(1:l_max,-l_max:l_max) C common /rexcf/ gelm,lpol C fac = dsqrt(x*x-1.d00) C do 200 l=0,ll_max,1 do 210 m=(-ll_max),ll_max,1 lpol(l,m) = 0.d00 210 continue 200 continue C lpol(0,0) = 1.d00 lpol(1,0) = x C do 100 l=2,ll_max,1 lpol(l,0) = (dfloat(2*l-1)*x*lpol(l-1,0) & -dfloat(l-1)*lpol(l-2,0))/dfloat(l) c write(6,*) l,lpol(l,0) 100 continue c write(6,*) C do 110 m=1,(ll_max-1),1 lpol(m,m) = (x*lpol(m,m-1) & -dfloat(2*m-1)*lpol(m-1,m-1))/fac lpol(m+1,m) = (2.d00*x*lpol(m+1,m-1) & -dfloat(2*m)*lpol(m,m-1))/fac do 120 l=(m+2),ll_max,1 lpol(l,m) = (dfloat(2*l-1)*x*lpol(l-1,m) & -dfloat(l+m-1)*lpol(l-2,m))/dfloat(l-m) 120 continue lpol(m,-m) = (x*lpol(m,1-m)-lpol(m-1,1-m)) & /(dfloat(2*m)*fac) lpol(m+1,-m) = (x*lpol(m+1,1-m)-lpol(m,1-m)) & /(dfloat(2*m+1)*fac) do 130 l=(m+2),ll_max,1 lpol(l,-m) = (dfloat(2*l-1)*x*lpol(l-1,-m) & -dfloat(l-m-1)*lpol(l-2,-m))/dfloat(l+m) 130 continue cc lpol(l,-1) = (x*lpol(l,0)-lpol(l-1,0))/(dfloat(l+1)*fac) c do 120 m=0,(l-1),1 c lpol(l,m+1) = (dfloat(l-m)*x*lpol(l,m) c & -dfloat(l+m)*lpol(l,m-1))/fac c 120 continue c do 130 m=(-1),(-l+1),(-1) c lpol(l,m-1) = (dfloat(l-m)*x*lpol(l,m) c & -fac*lpol(l,m+1))/dfloat(l+m) c 130 continue 110 continue m = ll_max lpol(m,m) = (x*lpol(m,m-1) & -dfloat(2*m-1)*lpol(m-1,m-1))/fac lpol(m,-m) = (x*lpol(m,1-m)-lpol(m-1,1-m)) & /(dfloat(2*m)*fac) C c write(6,*) 'legpol' c do 150 l=0,ll_max,1 c do 140 m=l,(-l),(-1) c write(6,*) l,m,lpol(l,m) c 140 continue c 150 continue C return end C ************************************* subroutine get_dfeldonr(xi,eps,lamax) C ************************************* implicit none integer*4 l_max,l,m,lamax parameter (l_max=3) real*8 xi,eps,fac,sum,res,ilm & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,ylmpih(0:l_max,-l_max:l_max),dfeldo(l_max) complex*16 i,selm(l_max,-l_max:l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /spheric/ ylmpih common /excf/ selm,dfeldo C do 100 l=1,l_max,1 fac = tpi*eps*eps/dfloat(2*l+1) fac = fac*fac/dfloat(2*l+1) res = 0.d00 do 110 m=l,(-l),(-2) C ++++++++++++++++++++++++++++++++++ call get_ilm(l,m,xi,eps,ilm,lamax) C ++++++++++++++++++++++++++++++++++ sum = ilm*ylmpih(l,m) res = res+sum*sum 110 continue dfeldo(l) = res*fac 100 continue C return end C ***************************************** subroutine get_dfeldor(e,v,a,theta,lamax) C ***************************************** implicit none integer*4 l_max,l,m,ll_max,lamax parameter (l_max=3,ll_max=l_max+1) real*8 e,v,a,theta,eps,gamma,b,k,xi,vi,fac,res,sum & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,dfeldo(l_max),kn(0:l_max) & ,lpol(0:ll_max,-ll_max:ll_max) complex*16 i,selm(l_max,-l_max:l_max) & ,gelm(1:l_max,-l_max:l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /excf/ selm,dfeldo common /bessel/ kn common /rexcf/ gelm,lpol C eps = 1.d00/dsin(theta*radfac/2.d00) gamma = 1.d00/dsqrt(1.d00-v*v) b = a*(dsqrt(eps*eps-1.d00)+pih/gamma) k = e/hbarc xi = k*b/(gamma*v) C C +++++++++++++++++++++ call get_besselkn(xi) C +++++++++++++++++++++ vi = 1.d00/v C +++++++++++++++++ call get_gelm(vi) C +++++++++++++++++ C fac = eps*eps/(2.d00*gamma) fac = fac*fac do 100 l=1,l_max,1 fac = fac*k*k*a*a res = 0.d00 if (l.gt.lamax) then dfeldo(l) = 0.d00 else do 110 m=l,(-l),(-1) sum = dreal(gelm(l,m)*dconjg(gelm(l,m))) & *kn(iabs(m))*kn(iabs(m)) res = res+sum 110 continue dfeldo(l) = res*fac endif 100 continue C return end C ************************** subroutine get_drot(theta) C ************************** implicit none integer*4 l_max,i1,i2,m,mp,j1,j2,l parameter (l_max=3) real*8 pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,theta,fac,cth,sth,clego,dl,dm,dmp,di1,di2,dj1,dj2 & ,rmat(3,3) complex*16 i,vu(-1:1,3) & ,dmpm(l_max,-l_max:l_max,-l_max:l_max) C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /rotmat/ dmpm C fac = 1.d00/dsqrt(2.d00) vu( 1,1) = -fac vu( 1,2) = -i*fac vu( 1,3) = 0.d00 vu( 0,1) = 0.d00 vu( 0,2) = 0.d00 vu( 0,3) = 1.d00 vu(-1,1) = fac vu(-1,2) = -i*fac vu(-1,3) = 0.d00 C do 100 i1=1,3,1 do 110 i2=1,3,1 rmat(i1,i2) = 0.d00 110 continue 100 continue C fac = theta*radfac/2.d00 cth = dcos(fac) sth = dsin(fac) C rmat(1,1) = -cth rmat(1,2) = -sth rmat(3,1) = sth rmat(3,2) = -cth rmat(2,3) = -1.d00 C do 200 mp=(-1),1,1 do 210 m=(-1),1,1 dmpm(1,mp,m) = 0.d00 do 220 i1=1,3,1 do 230 i2=1,3,1 dmpm(1,mp,m) = dmpm(1,mp,m) & +vu(m,i1)*rmat(i1,i2)*dconjg(vu(mp,i2)) 230 continue 220 continue 210 continue 200 continue C do 300 l=2,l_max,1 dl = dfloat(l) do 310 mp=-l,l,1 dmp = dfloat(mp) do 320 m=-l,l,1 dm = dfloat(m) dmpm(l,mp,m) = 0.d00 do 330 i1=-1,1,1 j1 = mp-i1 if (iabs(j1).le.(l-1)) then di1 = dfloat(i1) dj1 = dfloat(j1) do 340 i2=-1,1,1 j2 = m-i2 if (iabs(j2).le.(l-1)) then di2 = dfloat(i2) dj2 = dfloat(j2) dmpm(l,mp,m) = dmpm(l,mp,m) & +clego((dl-1.d00),dj1,1.d00,di1,dl,dmp) & *clego((dl-1.d00),dj2,1.d00,di2,dl,dm) & *dmpm(1,i1,i2)*dmpm(l-1,j1,j2) endif 340 continue endif 330 continue 320 continue 310 continue 300 continue C return end C *************************************** subroutine get_selmnr(xi,eps,v,a,lamax) C *************************************** implicit none integer*4 l_max,l,m,i1,mp,lamax parameter (l_max=3) real*8 xi,eps,v,a,ilm,fac & ,ylmpih(0:l_max,-l_max:l_max) & ,tmp(1:l_max,-l_max:l_max),dfeldo(1:l_max) complex*16 dmpm(l_max,-l_max:l_max,-l_max:l_max) & ,selm(l_max,-l_max:l_max) C common /spheric/ ylmpih common /excf/ selm,dfeldo common /rotmat/ dmpm C do 100 l=1,l_max,1 fac = v do 120 i1=1,l,1 fac = fac*a 120 continue do 110 m=l,(-l),(-2) C ++++++++++++++++++++++++++++++++++ call get_ilm(l,m,xi,eps,ilm,lamax) C ++++++++++++++++++++++++++++++++++ tmp(l,m) = ilm*ylmpih(l,m)/fac 110 continue do 130 m=(l-1),(-l+1),(-2) tmp(l,m) = 0.d00 130 continue 100 continue C do 200 l=1,l_max,1 do 210 m=-l_max,l_max,1 selm(l,m) = 0.d00 do 220 mp=-l_max,l_max,1 selm(l,m) = selm(l,m)+tmp(l,mp)*dmpm(l,mp,m) 220 continue 210 continue 200 continue C return end C ********************** subroutine init_ylmpih C ********************** implicit none integer*4 l_max,l,m,i1 parameter (l_max=3) real*8 fac,res,ylmpih(0:l_max,-l_max:l_max) & , pi,tpi,fpi,pih,radfac,e2,hbarc,amu complex*16 i C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu common /spheric/ ylmpih C do 100 l=0,l_max,1 do 110 m=(-l+1),(l-1),2 ylmpih(l,m) = 0.d00 110 continue fac = dfloat(2*l+1)/fpi do 120 m=(-l),l,2 res = 1.d00 do 130 i1=(l-m),1,(-1) res = res*dfloat(i1) 130 continue do 140 i1=(l+m),1,(-1) res = res*dfloat(i1) 140 continue res = dsqrt(res*fac) do 150 i1=(l-m),1,(-2) res = res/dfloat(i1) 150 continue do 160 i1=(l+m),1,(-2) res = res/dfloat(i1) 160 continue if (mod(((l+m)/2),2).ne.0) res=-res ylmpih(l,m) = res 120 continue 100 continue C return end C **************************************** subroutine get_ilm(l,m,xi,eps,ilm,lamax) C **************************************** implicit none integer*4 n_max parameter (n_max=5000) integer*4 l,m,lamax & ,ni1,ni2,i1,i2 real*8 xi,eps,ilm & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu & ,weps,xieps,dt,t,cht,sht,tmp(0:n_max),res,film,dum complex*16 i,tmp1,tmp2,tmp3 C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu C ilm = 0.d00 if (l.gt.lamax) return C weps = dsqrt(eps*eps-1.d00) xieps = xi*eps C ni1 = 40 C dt = 1.d07**(1.d00/dfloat(l))/eps dt = dlog(dt+dsqrt(dt*dt+1.d00)) if (xi.ne.0.d00) then dum = 1.d00+14.00/xieps dum = dlog(dum+dsqrt(dum*dum-1.d00)) if (dum.lt.dt) dt=dum endif ni2 = idnint(eps*dt) dt = 1.d00/(eps*dfloat(ni1)) C do 100 i1=0,ni1,1 t = dfloat(i1-ni1/2)*dt tmp(i1) = film(l,m,t,eps,weps,xi,xieps) 100 continue C ++++++++++++++++++++++++++++ call simpson(ni1,tmp,dt,res) C ++++++++++++++++++++++++++++ ilm = res do 110 i2=1,(ni2-1),1 do 120 i1=0,ni1,1 t = dfloat(i1+(2*i2-1)*(ni1/2))*dt tmp(i1) = film(l,m,t,eps,weps,xi,xieps) 120 continue C ++++++++++++++++++++++++++++ call simpson(ni1,tmp,dt,res) C ++++++++++++++++++++++++++++ ilm = ilm+2.d00*res 110 continue C ilm = dexp(-xi*pih)*ilm C return end C ******************************************************* double precision function film(l,m,t,eps,weps,xi,xieps) C ******************************************************* implicit none integer*4 l,m,i1 real*8 t,eps,weps,xi,xieps,cht,sht & ,pi,tpi,fpi,pih,radfac,e2,hbarc,amu complex*16 i,tmp1,tmp2,tmp3 C common /const/ i,pi,tpi,fpi,pih,radfac,e2,hbarc,amu C cht = dcosh(t) sht = dsinh(t) tmp1 = eps-weps*cht+i*sht tmp2 = 1.d00+i*eps*sht tmp3 = cdexp(i*xi*t) do 110 i1=1,l,1 tmp3 = tmp3/tmp2 110 continue tmp1 = tmp1/tmp2 do 120 i1=1,m,1 tmp3 = tmp3*tmp1 120 continue do 130 i1=(-1),m,(-1) tmp3 = tmp3/tmp1 130 continue film = dexp(-xieps*cht)*dreal(tmp3) C return end C ******************************* subroutine simpson(n,tmp,h,res) C ******************************* implicit none integer*4 n_max parameter (n_max=5000) integer*4 n,nn,i1 real*8 h,res,res1,res2,tmp(0:n_max) C if (n.gt.1) then nn = n-mod(n,2) res1 = 0.d00 res2 = 0.d00 do 100 i1=0,nn,2 res1 = res1+tmp(i1) 100 continue do 110 i1=1,(nn-1),2 res2 = res2+tmp(i1) 110 continue C res = (2.d00*res1+4.d00*res2-tmp(0)-tmp(nn))*h/3.d00 if (mod(n,2).eq.1) then res = res+h*(tmp(nn)+tmp(nn+1))/2.d00 endif else if (n.eq.1) then res = h*(tmp(0)+tmp(1))/2.d00 else res = 0.d00 endif endif C return end C ************************** subroutine get_besselkn(x) C ************************** implicit none integer*4 l_max,n parameter (l_max=3) real*8 x,k0,k1,kn(0:l_max) C common /bessel/ kn C C +++++++++++++++++++++++ call get_besselk0(x,k0) C +++++++++++++++++++++++ kn(0) = k0 C C +++++++++++++++++++++++ call get_besselk1(x,k1) C +++++++++++++++++++++++ kn(1) = k1 C do 100 n=2,l_max,1 kn(n) = kn(n-1)*dfloat(2*n-2)/x+kn(n-2) 100 continue C c do 200 n=0,l_max,1 c write(6,*) n,dexp(x)*kn(n) c 200 continue C return end C ***************************** subroutine get_besseli0(x,i0) C ***************************** implicit none integer*4 n real*8 x,i0,arg,sum,pi0s(0:6),pi0l(0:8) C data pi0s / 1.d00 & , 3.5156229d00, 3.0899424d00, 1.2067492d00 & , 0.2659732d00, 0.0360768d00, 0.0045813d00 / data pi0l / 0.39894228d00 & , 0.01328592d00, 0.00225319d00, -0.00157565d00 & , 0.00916281d00, -0.02057706d00, 0.02635537d00 & ,-0.01647633d00, 0.00392377d00 / C i0 = 0.d00 if (x.lt.0.d00) return if (x.lt.3.75d00) then arg = x*x/14.0625d00 sum = 0.d00 do 100 n=6,0,(-1) sum = arg*sum + pi0s(n) 100 continue i0 = sum else arg = 3.75d00/x sum = 0.d00 do 200 n=8,0,(-1) sum = arg*sum + pi0l(n) 200 continue i0 = sum*dexp(x)/dsqrt(x) endif C return end C ***************************** subroutine get_besseli1(x,i1) C ***************************** implicit none integer*4 n real*8 x,i1,arg,sum,pi1s(0:6),pi1l(0:8) C data pi1s / 0.5d00 & , 0.87890594d00, 0.51498869d00, 0.15084934d00 & , 0.02658733d00, 0.00301532d00, 0.00032411d00 / data pi1l / 0.39894228d00 & ,-0.03988024d00, -0.00362018d00, 0.00163801d00 & ,-0.01031555d00, 0.02282967d00, -0.02895312d00 & , 0.01787654d00, -0.00420059d00 / C i1 = 0.d00 if (x.lt.0.d00) return if (x.lt.3.75d00) then arg = x*x/14.0625d00 sum = 0.d00 do 100 n=6,0,(-1) sum = arg*sum + pi1s(n) 100 continue i1 = sum*x else arg = 3.75d00/x sum = 0.d00 do 200 n=8,0,(-1) sum = arg*sum + pi1l(n) 200 continue i1 = sum*dexp(x)/dsqrt(x) endif C return end C ***************************** subroutine get_besselk0(x,k0) C ***************************** implicit none integer*4 n real*8 x,k0,i0,arg,sum,pk0s(0:6),pk0l(0:6) C data pk0s / -0.57721566d00 & , 0.42278420d00, 0.23069756d00, 0.03488590d00 & , 0.00262698d00, 0.00010750d00, 0.00000740d00 / data pk0l / 1.25331414d00 & ,-0.07832358d00, 0.02189568d00, -0.01062446d00 & , 0.00587872d00, -0.00251540d00, 0.00053208d00 / C k0 = 0.d00 if (x.le.0.d00) return if (x.lt.2.d00) then arg = x*x/4.d00 sum = 0.d00 do 100 n=6,0,(-1) sum = arg*sum + pk0s(n) 100 continue C +++++++++++++++++++++++ call get_besseli0(x,i0) C +++++++++++++++++++++++ k0 = sum-dlog(x/2.d00)*i0 else arg = 2.d00/x sum = 0.d00 do 200 n=6,0,(-1) sum = arg*sum + pk0l(n) 200 continue k0 = sum*dexp(-x)/dsqrt(x) endif C return end C ***************************** subroutine get_besselk1(x,k1) C ***************************** implicit none integer*4 n real*8 x,k1,i1,arg,sum,pk1s(0:6),pk1l(0:6) C data pk1s / 1.d00 & , 0.15443144d00, -0.67278579d00, -0.18156897d00 & ,-0.01919402d00, -0.00110404d00, -0.00004686d00 / data pk1l / 1.25331414d00 & , 0.23498619d00, -0.03655620d00, 0.01504268d00 & ,-0.00780353d00, 0.00325614d00, -0.00068245d00 / C k1 = 0.d00 if (x.le.0.d00) return if (x.lt.2.d00) then arg = x*x/4.d00 sum = 0.d00 do 100 n=6,0,(-1) sum = arg*sum + pk1s(n) 100 continue C +++++++++++++++++++++++ call get_besseli1(x,i1) C +++++++++++++++++++++++ k1 = sum/x+dlog(x/2.d00)*i1 else arg = 2.d00/x sum = 0.d00 do 200 n=6,0,(-1) sum = arg*sum + pk1l(n) 200 continue k1 = sum*dexp(-x)/dsqrt(x) endif C return end C ------------------------------------------------ SUBROUTINE DFCOUL(L,ETA,RO,F,FP,G,GP,IEXP,SIGMA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION F(1),FP(1),G(1),GP(1),IEXP(1),SIGMA(1) DATA DEPI/6.283185307179586D0/ ETAC=ETA*ETA L1=L+1 CALL DFCZ0(ETA,RO,F0,FP0,G0,GP0,I,S) F(1)=F0 FP(1)=FP0 G(1)=G0 GP(1)=GP0 IEXP(1)=I SIGMA(1)=S IF(L)1,1,2 1 RETURN 2 LINF=0 IND=0 IF((ETA.GT.0).AND.(RO.LT.(ETA+ETA)))GO TO 21 Z=ETA+DSQRT(ETAC+DBLE(L*(L+1))) IF(RO.GE.Z)GO TO 20 7 ROINF=ETA+DSQRT(ETAC+DBLE(LINF*(LINF+1))) IF(RO-ROINF)3,4,4 4 IF(LINF-L)5,6,6 5 LINF=LINF+1 GO TO 7 3 IND=1 6 LIN=LINF+1 20 XM=1.D0 IF(IND.EQ.0)LIN=L1 DO 8 J=2,LIN ZIG=(DSQRT(ETAC+XM*XM))/XM ZAG=ETA/XM+XM/RO F(J)=(ZAG*F(J-1)-FP(J-1))/ZIG FP(J)=ZIG*F(J-1)-ZAG*F(J) G(J)=(ZAG*G(J-1)-GP(J-1))/ZIG GP(J)=ZIG*G(J-1)-ZAG*G(J) IEXP(J)=I SIG=SIGMA(J-1)+DATAN(ETA/(J-1)) IPI=SIG/DEPI SIG=SIG-IPI*DEPI IF(SIG)60,50,70 60 IF(SIG.LT.(-DEPI/2.D0))SIG=SIG+DEPI GO TO 50 70 IF(SIG.GT.(DEPI/2.D0))SIG=SIG-DEPI 50 SIGMA(J)=SIG 8 XM=XM+1.D0 IF(IND.EQ.0)RETURN GO TO 22 21 LIN=1 22 FTEST=F(LIN) FPTEST=FP(LIN) LMAX=LINF+25+IDINT(5.D0*DABS(ETA)) IF(LMAX-L)9,10,10 9 LMAX=L 10 FI=1.D0 FPI=1.D0 13 XM=LMAX+1 ZIG=(DSQRT(ETAC+XM*XM))/XM ZAG=ETA/XM+XM/RO FL=(ZAG*FI+FPI)/ZIG FPL=ZAG*FL-ZIG*FI IF(DABS(FL)-1.D15)26,27,27 26 IF(DABS(FPL)-1.D15)28,27,27 27 FL=FL*1.D-15 FPL=FPL*1.D-15 28 FI=FL FPI=FPL IF(LMAX-L)11,11,12 12 LMAX=LMAX-1 GO TO 13 11 F(LMAX+1)=FL FP(LMAX+1)=FPL IF(LMAX-LINF)15,15,14 14 GO TO 12 15 FACT=FTEST/F(LIN) FACTP=FPTEST/FP(LIN) INDICE=I/60 XM=LINF DO 16 J=LIN,L1 F(J)=F(J)*FACT FP(J)=FP(J)*FACTP 25 IF(J.EQ.1)GO TO 16 ZIG=(DSQRT(ETAC+XM*XM))/XM ZAG=ETA/XM+XM/RO G(J)=(ZAG*G(J-1)-GP(J-1))/ZIG GP(J)=ZIG*G(J-1)-ZAG*G(J) IF(DABS(G(J))-1.D60)17,18,18 17 IF(DABS(GP(J))-1.D60)19,18,18 18 G(J)=G(J)/1.D60 GP(J)=GP(J)/1.D60 INDICE=INDICE+1 19 IEXP(J)=INDICE*60 A=FP(J)*G(J) B=-F(J)*GP(J) IF(A-1.D0)29,30,30 29 I1=IDINT(DLOG10(A)) I2=IDINT(DLOG10(B)) GO TO 31 30 I1=IDINT(DLOG10(A))+1 I2=IDINT(DLOG10(B))+1 31 F(J)=F(J)*1.D1**(-I2) FP(J)=FP(J)*1.D1**(-I1) SIG=SIGMA(J-1)+DATAN(ETA/(J-1)) IPI=SIG/DEPI SIG=SIG-IPI*DEPI IF(SIG)61,51,71 61 IF(SIG.LT.(-DEPI/2.D0))SIG=SIG+DEPI GO TO 51 71 IF(SIG.GT.(DEPI/2.D0))SIG=SIG-DEPI 51 SIGMA(J)=SIG 16 XM=XM+1.D0 RETURN END SUBROUTINE DFCZ0(ETA,RO,F0,FP0,G0,GP0,IEXP,SIGMA0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION A1(110),A2(110),B1(110),B2(110) IF(RO)2,2,1 2 WRITE (6,1000) 1000 FORMAT (21H RO NEGATIF OU NUL **) RETURN 1 IF(ETA-30.D0)3,3,4 4 IF(DABS(ETA)-5.D2)28,28,29 28 CALL YFRICA(ETA,RO,F0,FP0,G0,GP0,IEXP,SIGMA0) 29 WRITE (6,1001) RETURN 1001 FORMAT (42H VALEUR ABSOLUE DE ETA SUPERIEURE A 500 **) 3 IF(ETA+8.D0)4,5,5 5 IF(ETA)6,7,6 7 F0=DSIN(RO) G0=DCOS(RO) FP0=G0 GP0=-F0 IEXP=0 SIGMA0=0.D0 RETURN 6 BORNE=1.666666666666667D0*DABS(ETA)+7.5D0 IF(RO-BORNE)8,9,9 9 CALL YFASYM(ETA,RO,F0,FP0,G0,GP0,IEXP,SIGMA0) RETURN 8 IF(ETA-10.D0)10,11,11 10 IF(ETA)12,12,13 13 IF(RO-2.D0)14,12,12 11 IF(ETA-(5.D0*RO+6.D1)/7.D0)12,12,14 12 CALL YFASYM(ETA,BORNE,F0,FP0,G0,GP0,IEXP,SIGMA0) H=BORNE DH=F0/H IF(ETA)20,21,21 20 N=-0.5D0*ETA+5.D0 GO TO 22 21 N=ETA/5.D0+5.D0 22 N=5*(N+1) Z=4.D0/H Y=1.D0-(ETA+ETA)*Z A1(N+2)=1.D-55 A1(N+3)=0.D0 A1(N+4)=1.D-64 B1(N+3)=1.D-50 B1(N+4)=1.D-68 A2(N+2)=0.D0 A2(N+3)=1.D-74 A2(N+4)=1.D-53 B2(N+3)=0.D0 B2(N+4)=1.D-66 M=N+2 DI=N DO 23 II=2,M I=M-II+2 B1(I)=B1(I+2)+ Z*(DI+1.D0)*A1(I+1) S=A1(I+2)+Y*(A1(I+1)-A1(I)) Q=(DI+2.D0)*B1(I)+(DI-1.D0)*B1(I+1) A1(I-1)=S-Z*Q B2(I)=B2(I+2)+ Z*(DI+1.D0)*A2(I+1) S=A2(I+2)+Y*(A2(I+1)-A2(I)) Q=(DI+2.D0)*B2(I)+(DI-1.D0)*B2(I+1) A2(I-1)=S-Z*Q IF(I.GE.N)GO TO 23 D=-(B2(I+2)+B2(I))/(B1(I+2)+B1(I)) DO 24 J=I,M A2(J)=A2(J)+D*A1(J) B2(J)=B2(J)+D*B1(J) 24 CONTINUE A2(I-1)=A2(I-1)+D*A1(I-1) 23 DI=DI-1.D0 Q=A1(3)-A1(1) C=A2(3)-A2(1) C=Q/C X1=0.5D0*(A1(2)-C*A2(2)) DO 25 I=3,M X1=X1+A1(I)-C*A2(I) 25 CONTINUE X1=DH/X1 X2=-C*X1 DO 26 I=2,M B1(I)=X1*B1(I)+X2*B2(I) A1(I)=X1*A1(I)+X2*A2(I) 26 CONTINUE A1(1)=X1*A1(1)+X2*A2(1) B1(1)=0.D0 X=RO/H Y=2.D0*X-1.D0 T1=1.D0 T2=Y RESULT=0.5D0*A1(2)+Y*A1(3) DERIVE=0.5D0*B1(2)+Y*B1(3) DO 27 I=2,N TI=2.D0*Y*T2-T1 T1=T2 T2=TI RESULT=RESULT+TI*A1(I+2) DERIVE=DERIVE+TI*B1(I+2) 27 CONTINUE F0=RESULT*RO FP0=DERIVE*RO+RESULT GO TO 30 C SERIE ORIGINE REGULIERE 14 PI=3.141592653589793D0 CALL JFLGAM(1.D0,ETA,TRA,SIGMA0,NTRUC) IEXP=0 RO2=RO*RO ETAP=ETA+ETA PIETA=PI*ETA PIETA2=0.5D0*PIETA B=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA) U0=0.D0 U1=RO U=U0+U1 UP=1.D0 XN=2.D0 DO 15 N=2,10000 XN1=XN*(XN-1.D0) U2=(ETAP*RO*U1-RO2*U0)/XN1 U=U+U2 UP=UP+XN*U2/RO 17 IF(DABS(U2/U).LT.1.D-10)GOTO19 18 U0=U1 U1=U2 XN=XN+1.D0 15 CONTINUE 19 F0=U/B FP0=UP/B 30 CALL YFIREG(ETA,RO,G0,GP0) RETURN END SUBROUTINE JFDELG (XD,YD,PAR,PAI,NBCHIF) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) DOUBLE PRECISION XD,YD,PAR,PAI,TEST,C,PI DOUBLE PRECISION X,Y,U,V,TRA,TRA1,COSI,SINI DOUBLE PRECISION COS2I,SIN2I,ZMOD,DEPI DOUBLE PRECISION TRB,XX DOUBLE PRECISION RAC2,PIS4 DOUBLE PRECISION SUPINT DIMENSION TEST(7),C(6) DATA TEST/2.9152D+7,2.2958D+3,1.4124D+2,3.9522D+1,19.6611D0,12.791 1D0,-10.D0/ DATA RAC2/0.3465735902799726/,PIS4/0.785398163397448/ DATA C/8.333333333333333D-2,-8.33333333333333D-3, 1 3.968253968253968D-3,-4.166666666666667D-3, 2 7.575757575757576D-3, -2.109279609279609D-2/ DATA PI/3.141592653589793/ DATA DEPI/6.283185307179586/ DATA SUPINT/2147483647.D0/ X=DABS(XD) XX=X NBCHIF=15 IF(YD)1,2,1 1 Y=DABS(YD) KR=1 I=DMOD(10.99D0-X,SUPINT) C TRANSLATION IF(I)3,3,4 4 TRA=I X=X+TRA C LOGARITHME(X+IY) (X,Y POSITIFS) 3 IF(X-Y)5,6,7 5 TRA1=X/Y TRB=1.D0+TRA1*TRA1 TRA=Y*DSQRT(TRB) SINI=1./(TRB*Y) COSI=SINI*TRA1 TRA1=Y/X GO TO 11 6 U=RAC2+DLOG(X) V=PIS4 SINI=0.5D0/X COSI=SINI GO TO 10 7 TRA1=Y/X TRB=1.D0+TRA1*TRA1 TRA=X*DSQRT(TRB) COSI=1./(TRB*X) SINI=COSI*TRA1 11 U=DLOG(TRA) V=DATAN(TRA1) C DEVELOPPEMENT ASYMPTOTIQUE ( X SUPERIEUR A 10 ) 10 PAR=U-0.5*COSI PAI=V+0.5*SINI ZMOD=X+Y IF(ZMOD-TEST(1))13,13,14 13 SIN2I=(SINI*COSI)+(SINI*COSI) COS2I=(COSI+SINI)*(COSI-SINI) SINI=SIN2I COSI=COS2I K=1 GO TO 15 16 TRA=COSI*COS2I-SINI*SIN2I SINI=SINI*COS2I+COSI*SIN2I COSI=TRA 15 PAR=PAR-C(K)*COSI PAI=PAI+C(K)*SINI K=K+1 IF(ZMOD-TEST(K))16,16,14 C TRANSLATION INVERSE 17 I=I-1 X=I X=XX+X 56 IF(X-Y)55,55,57 55 TRA1=X/Y TRB=X*TRA1+Y SINI=1.D0/TRB COSI=TRA1/TRB GO TO 19 57 TRA1=Y/X TRB=X+Y*TRA1 COSI=1.D0/TRB SINI=TRA1/TRB 19 PAR=PAR-COSI PAI=PAI+SINI 14 IF(I)18,18,17 C CONTROLE DU QUADRANT 18 IF(XD)20,61,21 61 TRA=PI*Y IF(TRA-1.D-2)300,300,301 300 TRB= TRA*(2.D0+TRA*(-2.D0+TRA*(1.333333333333333D0+ 1 TRA*( -0.6666666666666666D0+TRA*( 0.2666666666666666D0+ 2 TRA*( -0.08888888888888888D0+TRA* 0.02539682539682540D0 )))))) TRB=(2.D0-TRB)/TRB GO TO 302 301 TRB= DEXP(-TRA-TRA) TRB=(1.D0+TRB)/(1.D0-TRB) 302 PAI=0.5D0*(1.D0/Y+PI*TRB) 21 IF(YD)28,100,100 C X+IY CHANGE EN -X-IY 20 TRA=DEXP(-DEPI*Y) TRB=TRA*TRA COS2I=DEPI*DMOD(X,1.D0) SIN2I=-2.D0*TRA*DCOS(COS2I)+1.D0+TRB PAR=PAR+COSI+DEPI*TRA*DSIN(COS2I)/SIN2I PAI=PAI-SINI+PI*(TRB-1.D0)/SIN2I IF(YD)100,100,28 28 PAI=-PAI C ARGUMENT DANS -PI,PI 100 TRA=DABS(PAI/DEPI) IF(TRA-1.D+15)203,204,204 204 NBCHIF=0 PAI=0.D0 GO TO 29 203 IF(TRA-1.D0)205,205,206 206 NBCHIF=DLOG10(TRA) NBCHIF=14-NBCHIF TRA=DMOD(TRA,SUPINT) PAI=DMOD(TRA,1.D0)*DSIGN(DEPI,PAI) 205 IF(DABS(PAI)-PI)29,29,207 207 PAI=PAI-DSIGN(DEPI,PAI) GO TO 29 C DELGAMMA REEL 2 PAI=0.D0 IF(XD)31,32,33 C CONDITIONS D EXISTENCE 32 WRITE (6,1000) 1000 FORMAT (21H JFDELG(0) EST INFINI) GO TO 50 31 IF(X-4503599627370496.D0)34,35,35 35 WRITE (6,1001) 1001 FORMAT (30H ARGUMENT DE JFDELG TROP GRAND) GO TO 50 34 Y=DMOD(X,SUPINT) IF(Y) 400,36,400 400 IF(Y-0.99D0) 33,33,405 405 TRA= IDINT(Y+0.1D0) IF(DABS(Y-TRA)-5.D-15)36,36,33 36 WRITE (6,1002) 1002 FORMAT (28H JFDELG (-ENTIER) EST INFINI) 50 PAR=1.D+74 NBCHIF=0 GO TO 29 C TRANSLATION 33 I=DMOD(10.99D0-X,SUPINT) IF(I)37,37,38 38 TRA=I X=X+TRA C DEVELOPPEMENT ASYMPTOTIQUE 37 Y=DLOG(X) PAR=Y-0.5D0/X IF(X-TEST(1))39,39,43 39 COS2I=1.D0/(X*X) COSI=COS2I K=1 GO TO 41 42 COSI=COSI*COS2I 41 PAR=PAR-C(K)*COSI K=K+1 IF(X-TEST(K))42,42,40 C TRANSLATION INVERSE 44 I=I-1 X=I X=XX+X PAR=PAR-1.D0/X 40 IF(I)43,43,44 C X NEGATIF 43 IF(XD)45,29,29 45 PAR=PAR+1.D0/X Y=PI*DMOD(X,2.D0) PAR=PAR+PI*DCOS(Y)/DSIN(Y) 29 RETURN END SUBROUTINE YFASYM(ETA,RAU,FO,FPO,GO,GPO,IEXP,SIGO) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) IEXP=0 TRB=0.D0 RAU2=RAU+RAU ETAC=ETA*ETA CALL JFLGAM(1.D0,ETA,TRA,SIGO,NTRUC) 40 N=0 PS=1.D0 GS=0.D0 PT=0.D0 GT=1.D0-ETA/RAU SF=PS SG=GS SPF=PT SPG=GT 45 DENOM= DBLE(N+1)*RAU2 AN= DBLE(N+N+1)*ETA/DENOM BN= (ETAC- DBLE( N*(N+1)))/DENOM PS1=AN*PS-BN*PT GS1=AN*GS-BN*GT-PS1/RAU PT1=AN*PT+BN*PS GT1=AN*GT+BN*GS-PT1/RAU 42 SF=SF+PS1 SG=SG+GS1 SPF=SPF+PT1 SPG=SPG+GT1 N=N+1 IF(N-17)46,48,44 48 TRA=PS*PS+PT*PT 44 TRB=PS1*PS1+PT1*PT1 TEST=TRA-TRB IF(TEST)47,47,46 46 PS=PS1 GS=GS1 PT=PT1 GT=GT1 TRA=TRB GOTO 45 47 TETAO= RAU-ETA*DLOG (RAU2)+SIGO TRA= DSIN(TETAO) TRB=DCOS(TETAO) GO=SF*TRB-SPF*TRA GPO=SG*TRB-SPG*TRA FO=SPF*TRB+SF*TRA FPO=SPG*TRB+SG*TRA RETURN END SUBROUTINE YFIREG(ETA,RO,G0,GP0) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) IF(ETA.LE.0.D0)GOTO250 IF(ETA.LE.3.D0)GOTO251 IF(ETA.LE.1.D1)GOTO252 IF(ETA.LE.18.D0)GOTO253 IF(ETA.LE.22.D0)GOTO254 IF(RO.LE.0.3D0+(3.D1-ETA)/8.D1)GOTO200 C SERIE DE TAYLOR DEPART RAU0 300 CONTINUE RAU0=1.666666666666667D0*DABS(ETA)+7.5D0 CALL YFASYM(ETA,RAU0,F0,FP0,G0,GP0,IEXP,SIGMA0) X=RAU0-RO X2=X*X X3=X*X2 UNR=1.D0/RAU0 ETR0=1.D0-2.D0*ETA*UNR U0=G0 U1=-X*GP0 U2=-0.5D0*ETR0*X2*U0 S=U0+U1+U2 V1=U1/X V2=2.D0*U2/X T=V1+V2 XN=3.D0 DO10N=3,10000 C N=N XN1=XN-1.D0 XN1=XN*XN1 U3=X*U2*UNR*(1.D0-2.D0/XN)-ETR0*U1*X2/XN1+X3*U0*UNR/XN1 S=S+U3 V3=XN*U3/X T=T+V3 16 IF(DABS(U3/S).GT.1.D-11)GO TO 11 IF(DABS(V3/T).LE.1.D-11)GO TO 20 11 U0=U1 U1=U2 U2=U3 XN=XN+1.D0 10 CONTINUE 20 G0=S GP0=-T RETURN C SERIE ORIGINE 200 CONTINUE PI=3.141592653589793D0 GA=0.577215664901533D0 ETA2=ETA*ETA RO2=RO*RO ETAP=ETA+ETA PIETA=PI*ETA PIETA2=0.5D0*PIETA B=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA) CALL JFDELG(1.D0,ETA,PAR,PAI,NB) C1=ETAP*(GA+GA+DLOG(2.D0)-1.D0+PAR) U0=0.D0 U1=RO V0=1.D0 V1=C1*RO U=U0+U1 V=V0+V1 UP=1.D0 VP=C1 XN=2.D0 DO 104N=2,10000 XN1=XN*(XN-1.D0) U2=(ETAP*RO*U1-RO2*U0)/XN1 U=U+U2 V2=(ETAP*RO*V1-RO2*V0-ETAP*(XN+XN-1.D0)*U2)/XN1 V=V+V2 UP=UP+XN*U2/RO VP=VP+XN*V2/RO 17 IF(DABS(U2/U).GT.1.D-14)GOTO18 IF(DABS(V2/V).LE.1.D-14)GOTO19 18 U0=U1 U1=U2 V0=V1 V1=V2 XN=XN+1.D0 104 CONTINUE 19 GP=V+ETAP*U*DLOG(RO) G0=B*GP GP0=B*(VP+ETAP*(UP*DLOG(RO)+U/RO)) RETURN 250 IF(RO.LE.0.5D0*ETA+9.D0)GOTO200 GOTO 300 251 IF(RO.LE.2.25D0+7.35D0*(3.D0-ETA))GOTO200 GOTO 300 252 IF(RO.LE.1.2D0+1.5D-1*(1.D1-ETA))GOTO200 GOTO 300 253 IF(RO.LE.0.6D0+0.75D-1*(18.D0-ETA))GOTO200 GOTO 300 254 IF(RO.LE.0.4D0+0.5D-1*(22.D0-ETA))GOTO200 GOTO 300 END SUBROUTINE YFCLEN(ETA,RO,U,UP,V,VP,SIGMA0,IDIV,NN) IMPLICIT COMPLEX*16(A-D,F-H), REAL*8(E,P-Z) IMPLICIT INTEGER*4(I-N) C IF(NN.EQ.1)GO TO 20 C ETA2=ETA*ETA FA=DCMPLX(1.D0,ETA) M=0.25*ETA+4.D1 C C POLYNOMES DE TCHEBICHEV JUSQU'AU RANG M C K=M+2 X=(ETA+ETA)/RO XX=X+X-1.D0 T0=1.D0 T1=XX XX=XX+XX DO 6 J=2,K TJ=XX*T1-T0 T0=T1 T1=TJ 6 CONTINUE TM=T1 TL=T0 C C INITIALISATION C AM=(0.D0,0.D0) AL=(1.D-40,1.D-40) BN=(0.D0,0.D0) BM=(1.D-40,1.D-40) BL=(0.D0,0.D0) BK=DCMPLX(4.D0*DBLE(M+1),0.D0)*AL+BM F=(0.D0,0.D0) G =(0.D0,0.D0) GP=(0.D0,0.D0) C C RECURRENCE DESCENDANTE C K=M 10 R=K TK=XX*TL-TM TM=TL TL=TK HK=DCMPLX(TK,0.D0) C1=DCMPLX(R*(R+1.D0)-ETA2,ETA*(R+R+1.D0)) C2=(4.D0,0.D0)*DCMPLX(R+1.D0,0.D0) C2=C2*DCMPLX(-R-1.D0,ETA*3.D0) C3=FA*DCMPLX(-R-R-4.D0,ETA) C4=DCMPLX((7.D0*R+5.D0)/4.D0,0.D0) C5=DCMPLX(R+R+2.D0,0.D0) C6=DCMPLX((R+3.D0)/4.D0,0.D0) AK=(C2*AL+C3*AM-C4*BL-C5*BM-C6*BN)/C1 J=K/2 J=K-J-J IF(J)1,2,1 1 F=F-AK GO TO 3 2 F=F+AK 3 Z=CDABS(AK) G=G+HK*AK GP=GP+HK*BK C C F=A0/2-A1+A2-A3+A4-A5+... C C CONGRUENCE MODULO 10**60 C IF(Z-1.D60)4,5,5 5 D=(1.D60,0.D0) AK=AK/D AL=AL/D AM=AM/D BK=BK/D BL=BL/D BM=BM/D BN=BN/D F=F/D G=G/D GP=GP/D 4 IF(K)8,8,9 9 D=DCMPLX(4.D0*R,0.D0) BJ=D*AK+BL AM=AL AL=AK BN=BM BM=BL BL=BK BK=BJ K=K-1 GO TO 10 C C NORMALISATION ET CALCUL DE Z(RO) C 8 D=(0.5D0,0.D0)*AK F=F-D G=G-D GP=GP-(0.5D0,0.D0)*BK D=DCMPLX(0.D0,-ETA*DLOG(2.D0)+SIGMA0) AXPO=CDEXP(D) F=F/AXPO G=G/F GP=GP/F C C CALCUL DE F ET G C D=DCMPLX(0.D0,RO-ETA*DLOG(RO)) AXPO=CDEXP(D) D=G*AXPO GP=AXPO*(DCMPLX(0.D0,1.D0-ETA/RO)*G-DCMPLX(X/RO,0.D0)*GP) V=D D=(0.D0,-1.D0)*D U=D VP=GP GP=(0.D0,-1.D0)*GP UP=GP IDIV=0 RETURN C C SERIE ORIGINE C 20 PI=3.141592653589793D0 XA=0.577215664901533D0 RO2=RO*RO ETAP=ETA+ETA PIETA=PI*ETA Z=138.15510557964276D0 IDIV=0 IF(DABS(PIETA)-Z)21,21,22 22 INDG=IDINT(PIETA/Z) IDIV=60*INDG IF(ETA.LT.0) IDIV=0 PIETA=PIETA-Z*DBLE(INDG) 21 PIETA2=0.5D0*PIETA P=DEXP(PIETA2)*DSQRT(DSINH(PIETA)/PIETA) CALL JFDELG(1.D0,ETA,PAR,PAI,NB) Z1=ETAP*(XA+XA+DLOG(2.D0)-1.D0+PAR) U0=0.D0 U1=RO V0=1.D0 V1=Z1*RO U=U0+U1 V=V0+V1 UP=1.D0 VP=Z1 XN=2.D0 DO 104N=2,10000 XN1=XN*(XN-1.D0) U2=(ETAP*RO*U1-RO2*U0)/XN1 U=U+U2 V2=(ETAP*RO*V1-RO2*V0-ETAP*(XN+XN-1.D0)*U2)/XN1 V=V+V2 UP=UP+XN*U2/RO VP=VP+XN*V2/RO IF(DABS(U2/U).GT.1.D-14)GOTO18 IF(DABS(V2/V).LE.1.D-14)GOTO19 18 U0=U1 U1=U2 V0=V1 V1=V2 XN=XN+1.D0 104 CONTINUE 19 PP=V+ETAP*U*DLOG(RO) W=U/P WP=UP/P V=P*PP VP=P*(VP+ETAP*(UP*DLOG(RO)+U/RO)) U=W UP=WP RETURN END SUBROUTINE YFRICA(ETA,RAU,FO,FPO,GO,GPO,IDIV,SIGMA0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION Q(5),QP(5) C C COEFFICIENTS RICCATI C DATA G61,G62,G63,G64,G65,G66,G67,G68,G69, 1 G610,G611/ 2 0.1159057617187498D-01, 0.3863525390624998D-01, 3 0.4660034179687498D-01, 0.4858398437499998D-01, 4 0.1156514485677080D 01, 0.5687475585937496D 01, 5 0.1323888288225445D 02, 0.1713083224826384D 02, 6 0.1269003295898436D 02, 0.5055236816406248D 01, 7 0.8425394694010415D 00/ DATA G81,G82,G83,G84,G85,G86,G87,G88,G89, 1 G810,G811,G812,G813,G814,G815/ 2 0.1851092066083633D-01, 0.8638429641723630D-01, 3 0.1564757823944092D 00, 0.1430139541625977D 00, 4 0.1924622058868408D 00, 0.8500803152720129D 01, 5 0.7265429720878595D 02, 0.3057942376817972D 03, 6 0.7699689544836672D 03, 0.1254157054424285D 04, 7 0.1361719536066055D 04, 0.9831831171035763D 03, 8 0.4547869927883148D 03, 0.1222640538215636D 03, 9 0.1455524450256709D 02/ DATA GP61,GP62,GP63,GP64,GP65,GP66/ 2 0.2897644042968748D-01, 0.2318115234375000D 00, 3 0.8056640625000000D 00, 0.1601562499999998D 01, 4 0.3046875000000000D 00, 0.5624999999999998D 01/ DATA GP81,GP82,GP83,GP84,GP85,GP86,GP87,GP88/ 1 0.6478822231292720D-01, 0.6910743713378906D 00, 2 0.3322952270507811D 01, 0.9483032226562498D 01, 3 0.1769653320312499D 02, 0.3478710937499998D 02, 4 0.5020312499999999D 02, 0.7874999999999999D 02/ DATA Q /0.4959570165D-1,0.8888888889D-2,0.2455199181D-2, $ 0.9108958061D-3,0.2534684115D-3/ DATA QP /0.1728260369D0,0.3174603174D-3,0.3581214850D-2, 1 0.3117824680D-3,0.9073966427D-3/ CALL JFLGAM(1.D0,ETA,TRA,SIGMA0,IND) RAU2=RAU+RAU RAUC=RAU*RAU ETAC=ETA*ETA ETA2=ETA+ETA ETARO=ETA*RAU ETARO2=ETARO+ETARO PIETA=3.141592653589793*ETA IND=0 JND=0 IG=0 IF(ETA)20,20,21 20 IF(-ETARO-14.0625D0)32,22,22 22 INDICE=1 C RICCATI 3 IDIV=0 GO TO 2 21 IF(DABS(RAU-ETA2).LE.1.D-9)GO TO 18 IF(RAU-ETA2)30,18,31 31 IF(RAU-ETA2-2.D1*(ETA**0.25D0))34,33,33 33 INDICE=0 C RICCATI 2 IDIV=0 GO TO 2 32 NN=1 GO TO 35 34 NN=0 35 CALL YFCLEN(ETA,RAU,FO,FPO,GO,GPO,SIGMA0,IDIV,NN) RETURN 30 IF(ETARO-12.D0)32,32,23 23 TRA=ETA2-6.75D0*(ETA**0.4D0) IF(RAU-TRA)6,6,24 24 IND=1 JND=1 RO=RAU RAU=TRA RAU0=TRA C RICCATI 1 6 X=RAU/ETA2 U= (1.D0-X)/X X2=X*X RU= DSQRT(U) RX= DSQRT(X) TRE= 1.D0/(U*RU*ETA2) TRB=TRE*TRE FI= (DSQRT((1.D0-X)*X)+ DASIN(RX)-1.570796326794897D0)*ETA2 TR1= -0.25D0*DLOG(U) 602 TR2= -((9.D0*U+6.D0)*U+5.D0)/48.D0 603 TR3= ((((-3.D0*U-4.D0)*U+6.D0)*U+12.D0)*U+5.D0)/64.D0 604 TR4=- ((((((U+2.D0)*945.D0*U+1395.D0)*U+12300.D0)*U+25191.D0)*U 1 +19890.D0)*U+5525.D0)/46080.D0 605 TR5= (((((((( -27.D0*U-72.D0)*U-68.D0)*U+360.D0)*U+2190.D0)*U 1 +4808.D0)*U+5148.D0)*U+2712.D0)*U+565.D0)/2048.D0 606 TR6=- (((((((((G61*U+G62)*U+G63)*U+G64)*U+G65)*U+G66)*U+G67)*U 1 +G68)*U+G69)*U+G610)*U+G611 607 TR7= ((((((((((((-81.*U-324.)*U-486.)*U-404.)*U+4509.)*U+52344.) 1 *U+233436.)*U+567864.)*U+838521.)*U+775884.)*U 2 +441450.)*U+141660.)*U+19675.) /6144. 608 TR8= ((((((((((((( G81*U+G82)*U+G83)*U+G84)*U+G85)*U+G86)*U 1 +G87)*U+G88)*U+G89)*U+G810)*U+G811)*U+G812)*U 2 +G813)*U+G814)*U+G815 FI= FI+TRE*(TR2+TRB*(TR4+TRB*(TR6+TRB*TR8))) PSI=-FI TRA= TR1+TRB*(TR3+TRB*(TR5+TRB*TR7)) FI=FI+TRA PSI=PSI+TRA FIP=RU*ETA2 TRA=1.D0/(X2*U) TR1=0.25D0 TRE= TRE/(X2*X2*U) TRB=TRB/(X2*X2) 702 TR2= -(8.D0*X-3.D0)/32.D0 703 TR3= (( 24.D0*X-12.D0)*X+3.D0)/64.D0 704 TR4= ((( -1536.D0*X+704.D0)*X-336.D0)*X+63.D0)/2048.D0 705 TR5= (((( 1920.D0*X-576.D0)*X+504.D0)*X-180.D0)*X+27.D0)/1024.D0 706 TR6 = (((( -GP66*X+GP65)*X-GP64)*X+GP63)*X-GP62)*X+GP61 707 TR7= - (((((( -40320.D0*X-10560.D0)*X-13248.D0)*X+7560.D0)*X 1 -3132.D0)*X+756.D0)*X-81.D0) / 2048.D0 708 TR8 =- (((((( GP88*X+GP87)*X+GP86)*X-GP85)*X+GP84)*X-GP83)*X+GP82 1 )*X-GP81 FIP=FIP +TRE*(TR2+TRB*(TR4+TRB*(TR6+TRB*TR8))) PSIP=-FIP TRA= TRA*(TR1+TRB*(TR3+TRB*(TR5+TRB*TR7))) FIP=FIP+TRA PSIP=PSIP+TRA XXX=138.1551055796428D0 INDG=IDINT(PSI/XXX) IDIV=60*INDG IF(INDG.EQ.0)GO TO 8 PSI=PSI-XXX*DBLE(INDG) FI =FI +XXX*DBLE(INDG) 8 FO=0.5D0*DEXP(FI) GO= DEXP( PSI) FPO= FO*FIP/ETA2 GPO=GO*PSIP/ETA2 IF(JND.EQ.0)RETURN RAU=RO GO=FO GPO=FPO 27 X=RAU0-RO X2=X*X X3=X*X2 UNR=1.D0/RAU0 ETR0=1.D0-2.D0*ETA*UNR U0=GO U1=-X*GPO U2=-0.5D0*ETR0*X2*U0 S=U0+U1+U2 V1=U1/X V2=2.D0*U2/X T=V1+V2 XN=3.D0 DO10N=3,10000 C N=N XN1=XN-1.D0 XN1=XN*XN1 U3=X*U2*UNR*(1.D0-2.D0/XN)-ETR0*U1*X2/XN1+X3*U0*UNR/XN1 S=S+U3 V3=XN*U3/X T=T+V3 16 IF(DABS(U3/S).GT.1.D-10)GO TO 11 IF(DABS(V3/T).LE.1.D-10)GO TO 17 11 U0=U1 U1=U2 U2=U3 XN=XN+1.D0 10 CONTINUE 17 IF(IG)25,26,25 25 GO=S GPO=-T FO=HO FPO=HPO RETURN 26 HO=S HPO=-T 18 ET0=ETA**(0.166666666666667) ETAD=ETAC*ETAC ET=ETA**(0.6666666666666667) ET1=ET*ET ET2=ET1*ET1 ET3=ET2*ET ET4=ETAD*ET ET5=ET4*ET FO=1.D0-Q(1)/ET1-Q(2)/ETAC-Q(3)/ET3-Q(4)/ETAD-Q(5)/ET5 GO=1.D0+Q(1)/ET1-Q(2)/ETAC+Q(3)/ET3-Q(4)/ETAD+Q(5)/ET5 FPO=1.D0+QP(1)/ET+QP(2)/ETAC+QP(3)/ET2+QP(4)/ETAD+QP(5)/ET4 GPO=1.D0-QP(1)/ET+QP(2)/ETAC-QP(3)/ET2+QP(4)/ETAD-QP(5)/ET4 FO=0.7063326373D0*ET0*FO GO=1.223404016D0*ET0*GO FPO=0.4086957323D0*FPO/ET0 GPO=-0.7078817734D0*GPO/ET0 IDIV=0 IF(IND.EQ.0)RETURN IG=1 RAU0=ETA2 GO TO 27 2 X=ETA2/RAU X2=X*X U=1.D0-X RU= DSQRT(U) U3=U*U*U TRD= 1.D0/(U3*ETA2*ETA2) TRC=X2*TRD TRE=1.D0/(U*RU*ETA2) FI= -0.25D0*DLOG(U) TRB=TRD/64.D0 TR3= (((3.D0*U-4.D0)*U-6.D0)*U+12.D0)*U-5.D0 501 TR5= (((((((( -27.D0*U+72.D0)*U-68.D0)*U-360.D0)*U+2190.D0)*U 1 -4808.D0)*U+5148.D0)*U-2712.D0)*U+565.D0)/32.D0 502 TR7= (((((((((((( 81.D0*U-324.D0)*U+486.D0)*U-404.D0)*U-4509.D0)*U 1 +52344.D0)*U-233436.D0)*U+567864.D0)*U-838521.D0)*U+775884.D0 2 )*U-441450.D0)*U+141660.D0)*U-19675.D0)/96.D0 FI= FI+TRB*(TR3+TRD*(TR5+TRD*TR7)) FIP=0.25D0/U TRB=3.D0*TRC/(64.D0*U) TR3= (X-4.D0)*X+8.D0 511 TR5= ((((9.D0*X-60.D0)*X+168.D0)*X-192.D0)*X+640.D0)/16.D0 512 TR7= ((((((-27.D0*X+252.D0)*X-1044.D0)*X+2520.D0)*X-4416.D0)*X 1 -3520.D0)*X-13440.D0)/32.D0 FIP =FIP+TRB*(TR3+TRC*(TR5+TRC*TR7)) TRA= DABS((RU-1.D0)/(RU+1.D0)) PSI= (0.5D0*DLOG(TRA)+RU/X)*ETA2+0.785398163397448D0 TR2= -((9.D0*U-6.D0)*U+5.D0)/48.D0 521 TR4= ((((((U-2.D0)*945.D0*U+1395.D0)*U-12300.D0)*U+25191.D0)*U 1 -19890.D0)*U+5525.D0)/46080.D0 522 TR6 = (((((((((-G61*U+G62)*U-G63)*U+G64)*U-G65)*U+G66)*U 1 -G67)*U+G68)*U-G69)*U+G610)*U-G611 523 TR8= ((((((((((((( G81*U-G82)*U+G83)*U-G84)*U+G85)*U-G86)*U 1 + G87)*U-G88)*U+G89)*U-G810)*U+G811)*U 2 -G812)*U+G813)*U-G814)*U+G815 PSI= PSI+TRE*(TR2+TRD*(TR4+TRD*(TR6+TRD*TR8))) PSIP = -RU*ETA2/X2 TRB=TRE*X/U TR2= (3.D0*X-8.D0)/32.D0 531 TR4= - (((63.D0*X-336.D0)*X+704.D0)*X-1536.D0)/2048.D0 532 TR6 = (((( GP61*X-GP62)*X+GP63)*X-GP64)*X+GP65)*X-GP66 533 TR8 = (((((( -GP81*X+GP82)*X-GP83)*X+GP84)*X-GP85)*X+GP86)*X 1 +GP87)*X+GP88 PSIP =PSIP+ TRB*(TR2+TRC*(TR4+TRC*(TR6+TRC*TR8))) TRA= DEXP(FI) FO= TRA*DSIN(PSI) GO= TRA*DCOS(PSI) IF(INDICE)535,536,535 535 TRA=FO FO=GO GO=-TRA 536 TRA=-ETA2/RAUC FPO=(FIP*FO+PSIP*GO)*TRA GPO=(FIP*GO-PSIP*FO)*TRA RETURN END SUBROUTINE JFLGAM(XD,YD,PAR,PAI,NBCHIF) IMPLICIT INTEGER*4(I-N) DOUBLE PRECISION XD,YD,PAR,PAI,TEST,C,HLO2PI,PI,PIS2,PIS4 DOUBLE PRECISION X,Y,U,V,TRA,TRA1,ALO2PI,RAC2,COSI,SINI DOUBLE PRECISION COS2I,SIN2I,ZMOD,DEPI DOUBLE PRECISION XX DOUBLE PRECISION ALOPI DOUBLE PRECISION SUPINT DIMENSION TEST(7),C(6) DATA TEST/2.9152D+7,2.2958D+3,1.4124D+2,3.9522D+1,19.6611D0,12.791 1D0,-10.D0/ DATA C/8.333333333333333D-2,-2.777777777777777D-3,7.936507936507 1937D-4,-5.952380952380952D-4,8.417508417508418D-4,-1.917526917526 2918D-3/ DATA HLO2PI/0.918938533204672/,PI/3.141592653589793/ DATA PIS2/1.570796326794897/,PIS4/0.785398163397448/ DATA ALO2PI/1.837877066409345/,RAC2/0.3465735902799726/ DATA DEPI/6.283185307179586/,ALOPI/1.1447298858494002/ DATA SUPINT/2147483647.D0/ NBCHIF=15 X=DABS(XD) XX=X IF(YD)1,2,1 1 Y=DABS(YD) KR=1 I=DMOD(10.99D0-X,SUPINT) C TRANSLATION IF(I)3,3,4 4 TRA=I X=X+TRA C LOGARITHME(X+IY) (X,Y POSITIFS) 3 IF(X-Y)5,6,7 5 TRA1=X/Y IF(TRA1)8,8,9 8 U=DLOG(Y) V=PIS2 GO TO 10 6 U=RAC2+DLOG(X) V=PIS4 GO TO 10 9 TRA=Y*DSQRT(1.D0+TRA1*TRA1) TRA1=Y/X 11 U=DLOG(TRA) V=DATAN(TRA1) 10 GO TO (12,19,23),KR 7 TRA1=Y/X TRA=X*DSQRT(1.D0+TRA1*TRA1) GO TO 11 12 KR=2 C DEVELOPPEMENT ASYMPTOTIQUE ( X SUPERIEUR A 10 ) TRA=X-0.5D0 PAI=V*TRA+Y*(U-1.D0) PAR=-X+HLO2PI+U*TRA-V*Y ZMOD=X+Y IF(ZMOD-TEST(1))13,13,14 13 TRA=X*X+Y*Y COSI=X/TRA SINI=Y/TRA SIN2I=(SINI*COSI)+(SINI*COSI) COS2I=(COSI+SINI)*(COSI-SINI) K=1 GO TO 15 16 TRA=COSI*COS2I-SINI*SIN2I SINI=SINI*COS2I+COSI*SIN2I COSI=TRA 15 PAR=PAR+C(K)*COSI PAI=PAI-C(K)*SINI K=K+1 IF(ZMOD-TEST(K))16,16,14 C TRANSLATION INVERSE 17 I=I-1 X=I X=XX+X GO TO 3 19 PAR=PAR-U PAI=PAI-V 14 IF(I-1)18,60,17 60 IF(XD)17,61,17 C CONTROLE DU QUADRANT 18 IF(XD)20,61,21 61 TRA=PI*Y IF(TRA-1.D-2)300,300,301 300 PAR= TRA*(2.D0+TRA*(-2.D0+TRA*(1.333333333333333D0+ 1 TRA*( -0.6666666666666666D0+TRA*( 0.2666666666666666D0+ 2 TRA*( -0.08888888888888888D0+TRA* 0.02539682539682540D0 )))))) TRA1=-DLOG(Y)-DLOG(PAR) GO TO 302 301 PAR=1.D0-DEXP(-TRA-TRA) TRA1=-DLOG(Y*PAR) 302 PAR=0.5D0*(ALO2PI-TRA+TRA1) PAI=PAI-PIS2 21 IF(YD)28,100,100 C X+IY CHANGE EN -X-IY 20 TRA=PI*Y PAR=ALO2PI-U-PAR-TRA PAI=PI-V-PAI TRA=DEXP(-TRA-TRA) X=PI*DMOD(X,2.D0) SINI=(1.D0-TRA)*DCOS(X) COSI=(1.D0+TRA)*DSIN(X) KR=3 X=DABS(COSI) Y=DABS(SINI) GO TO 3 23 IF(COSI)24,25,25 24 V=PI-DSIGN(V,SINI) GO TO 26 25 IF(SINI)27,26,26 27 V=-V 26 PAR=PAR-U PAI=PAI-V IF(YD)100,100,28 28 PAI=-PAI C ARGUMENT DANS -PI,PI 100 TRA=DABS(PAI/DEPI) IF(TRA-1.D+15)203,204,204 204 NBCHIF=0 PAI=0.D0 GO TO 29 203 IF(TRA-1.D0)205,205,206 206 NBCHIF=DLOG10(TRA) NBCHIF=14-NBCHIF TRA=DMOD(TRA,SUPINT) PAI=DMOD(TRA,1.D0)*DSIGN(DEPI,PAI) 205 IF(DABS(PAI)-PI)29,29,207 207 PAI=PAI-DSIGN(DEPI,PAI) GO TO 29 C JFLGAM REEL 2 PAI=0.D0 IF(XD)31,32,33 C CONDITIONS D EXISTENCE 32 WRITE (6,1000) 1000 FORMAT (21H JFLGAM(0) EST INFINI) GO TO 50 31 IF(X-4503599627370496.D0)34,35,35 35 WRITE (6,1001) 1001 FORMAT (30H ARGUMENT DE JFLGAM TROP GRAND) GO TO 50 34 Y=DMOD(X,SUPINT) IF(Y)400,36,400 400 IF(Y-0.99D0)33,33,405 405 TRA=IDINT(Y+0.1D0) IF(DABS(Y-TRA)-5.D-15)36,36,33 36 WRITE (6,1002) 1002 FORMAT (28H JFLGAM (-ENTIER) EST INFINI) 50 PAR=1.D+74 NBCHIF=0 GO TO 29 C TRANSLATION 33 I=DMOD(10.99D0-X,SUPINT) IF(I)37,37,38 38 TRA=I X=X+TRA C DEVELOPPEMENT ASYMPTOTIQUE 37 Y=DLOG(X) PAR=-X+HLO2PI+Y*(X-0.5D0) IF(X-TEST(1))39,39,43 39 COSI=1.D0/X COS2I=COSI*COSI K=1 GO TO 41 42 COSI=COSI*COS2I 41 PAR=PAR+C(K)*COSI K=K+1 IF(X-TEST(K))42,42,40 C TRANSLATION INVERSE 44 X=X-1.D0 48 Y=DLOG(X) PAR=PAR-Y I=I-1 40 IF(I-1)43,49,44 49 X=DABS(XD) GO TO 48 C X NEGATIF 43 IF(XD)45,29,29 45 PAR=ALOPI-PAR-Y Y=PI*DMOD(X,2.D0) Y=-DSIN(Y) IF(Y)46,36,47 46 Y=-Y PAI=PI 47 PAR=PAR-DLOG(Y) 29 RETURN END C ******************************************** double precision function sjsym(a,b,e,d,c,f) C ******************************************** C Calculation of 6j-symbol | a b c | C | d e f | implicit none integer*4 i,minz,maxz real*8 a,b,c,d,e,f, & fak1,fak2,facu,arg,sum,min,max,z,tria,signum external facu,tria,signum C sjsym=0.0 fak1=tria(a,b,e)*tria(c,d,e)*tria(a,c,f)*tria(b,d,f) C min=-1.0 arg=a+b+e if (arg.gt.min) min=arg arg=c+d+e if (arg.gt.min) min=arg arg=a+c+f if (arg.gt.min) min=arg arg=b+d+f if (arg.gt.min) min=arg max=a+b+c+d arg=a+d+e+f if (arg.lt.max) max=arg arg=b+c+e+f if (arg.lt.max) max=arg C if (min.gt.max) min=max sum=signum(min)*facu(min+1.0) sum=sum/facu(min-a-b-e) sum=sum/facu(min-c-d-e) sum=sum/facu(min-a-c-f) sum=sum/facu(min-b-d-f) sum=sum/facu(a+b+c+d-min) sum=sum/facu(a+d+e+f-min) sum=sum/facu(b+c+e+f-min) fak2=sum minz=idnint(min) maxz=idnint(max) do 10 i=(minz+1),maxz,1 C z=dflotj(i) C 2001/06/25 z=dfloat(i) sum=-sum*(z+1.0) sum=sum*(a+b+c+d-z+1.0)*(a+d+e+f-z+1.0)*(b+c+e+f-z+1.0) sum=sum/((z-a-b-e)*(z-c-d-e)*(z-a-c-f)*(z-b-d-f)) fak2=fak2+sum 10 continue sjsym=fak1*fak2 return end C *********************************** double precision function signum(x) C *********************************** implicit none integer*4 i,ix,val real*8 x C ix=iabs(idnint(x)) val=1 do 10 i=1,ix,1 val=-val 10 continue C signum=dflotj(val) C 2001/06/25 signum=dfloat(val) return end C ************************************* double precision function tria(a,b,c) C ************************************* C Calculation of triangle function \Delta(a,b,c) C implicit none real*8 a,b,c,val,facu external facu C tria=0.0 if ((abs(a-b).lt.(c+0.1)).and.((a+b+0.1).gt.c)) then val=facu(a+b-c) else return endif if ((abs(a-c).lt.(b+0.1)).and.((a+c+0.1).gt.b)) then val=val*facu(a+c-b) else return endif if ((abs(b-c).lt.(a+0.1)).and.((b+c+0.1).gt.a)) then val=val*facu(b+c-a) else return endif val=val/facu(a+b+c+1.0) tria=dsqrt(val) return end C ************************************************ double precision function clego(j1,m1,j2,m2,j,m) C ************************************************ C Calculation of Clebsch-Gordan coefficient ( j1, m1, j2, m2 | j, m) C implicit none integer*4 i,min,max,narg real*8 j1,j2,j,m1,m2,m,fak1,fak2,fak3,facu, & arg,sum,minz,maxz external facu C clego=0.0 narg=(2.0*(j1+j2+j)+0.1) if (mod(narg,2).ne.0) return if (abs(m1).gt.(j1+0.1)) return if (abs(m2).gt.(j2+0.1)) return if (abs(m).gt.(j+0.1)) return if (abs(m1+m2-m).gt.0.1) then return else arg=j1+j2-j if ((abs(j1-j2).lt.(j+0.1)).and.((j1+j2+0.1).gt.j)) then fak1=facu(arg) else return endif arg=j+j1-j2 if ((abs(j-j1).lt.(j2+0.1)).and.((j+j1+0.1).gt.j2)) then fak1=fak1*facu(arg) else return endif arg=j+j2-j1 if ((abs(j-j2).lt.(j1+0.1)).and.((j+j2+0.1).gt.j1)) then fak1=fak1*facu(arg) else return endif arg=j1+j2+j+1 fak1=fak1/facu(arg) fak1=fak1*(2*j+1) C fak2=facu(j1+m1)*facu(j1-m1) fak2=fak2*facu(j2+m2)*facu(j2-m2) fak2=fak2*facu(j+m)*facu(j-m) C arg=j1+j2-j maxz=arg arg=j1-m1 if (arg.lt.maxz) maxz=arg arg=j2+m2 if (arg.lt.maxz) maxz=arg minz=0.0 arg=-(j-j2+m1) if (arg.gt.minz) minz=arg arg=-(j-j1-m2) if (arg.gt.minz) minz=arg C if (minz.gt.maxz) minz=maxz sum=1.0/facu(minz) arg=j1+j2-j-minz sum=sum/facu(arg) arg=j1-m1-minz sum=sum/facu(arg) arg=j2+m2-minz sum=sum/facu(arg) arg=j-j2+m1+minz sum=sum/facu(arg) arg=j-j1-m2+minz sum=sum/facu(arg) max=(maxz+0.1) min=(minz+0.1) fak3=sum C write(2,*) 'min ',min,' max ',max do 20 i=(min+1),max,1 sum=sum*(j1+j2-j-i+1)*(j1-m1-i+1)*(j2+m2-i+1) sum=-sum/(i*(j-j2+m1+i)*(j-j1-m2+i)) fak3=fak3+sum 20 continue C if (min.lt.0) min=-min do 30 i=1,min,1 fak3=-fak3 30 continue C write(6,*) fak1,fak2,fak3 clego=dsqrt(fak1*fak2)*fak3 endif return end C ********************************* double precision function facu(x) C ********************************* C Calculation of factorial function x! C implicit none integer*4 i,n real*8 x,j C if (x.lt.0.1) then facu=1.0 else j=1.0 n=idnint(x) do 10 i=1,n,1 C j=j*dflotj(i) C 2001/06/25 j=j*dfloat(i) 10 continue facu=j endif return end C ************************************************** double precision function njsym(a,b,c,d,e,f,g,h,i) C ************************************************** C Calculation of 9j-symbol | a b c | C | d e f | C | g h i | C implicit none real*8 k,k1,k2,kmin,kmax,a,b,c,d,e,f,g,h,i,sum,prod, & fak1,fak2,fak3,sjsym,signum external sjsym,signum C sum=0.0 k1=dabs(d-h) k2=dabs(b-f) if (k1.lt.k2) then kmin=k1 else kmin=k2 endif k1=d+h k2=b+f if (k1.gt.k2) then kmax=k1 else kmax=k2 endif do 10 k=kmin,(kmax+0.001),1.0 fak1=sjsym(a,d,g,h,i,k) if (fak1.ne.0.0) then fak2=sjsym(b,e,h,d,k,f) if (fak2.ne.0.0) then fak3=sjsym(c,f,i,k,a,b) endif endif prod=fak1*fak2*fak3*(2.0*k+1.0)*signum(2.0*k) sum=sum+prod 10 continue njsym=sum return end C ****************************************** subroutine simul(nev,file_alm,file_ev,acc) C ****************************************** implicit none integer*4 nt_max,ne_max,l_max,tl_max,dim parameter (nt_max=90,ne_max=501,l_max=3,tl_max=2*l_max,dim=5) integer*4 nev,file_alm,file_ev,nt,ne,i1,i2,l,m,n,nc,nd,k1,k2 real*8 m_x,m_a,m_b,m_c,e_i,t_p,t_min,t_max,e_min,e_max & ,dt,de,vt(0:nt_max),ve(0:ne_max),fac(0:nt_max,0:ne_max) & ,x(dim),d(dim),wi,work(dim),pb(4),pc(4),px(4) & ,e_rel,theta_cm,phi_cm,theta_rel,phi_rel,acc & ,pi,tpi,radfac complex*16 alm(0:nt_max,0:ne_max,0:tl_max,-tl_max:tl_max) C common /simpar0/ pi,tpi,radfac common /simpar1/ m_x,m_a,m_b,m_c,e_i,t_p common /simpar2/ t_min,t_max,dt,e_min,e_max,de,nt,ne common /simpar3/ vt,ve common /simpar4/ alm common /simpar5/ fac C pi = 4.d00*datan(1.d00) tpi = 2.d00*pi radfac = pi/180.d00 C read(file_alm,*) m_x,m_a read(file_alm,*) m_b,m_c read(file_alm,*) e_i,t_p read(file_alm,*) t_min,t_max read(file_alm,*) e_min,e_max read(file_alm,*) nt,ne C dt = (t_max-t_min)/dfloat(nt) do 100 i1=0,nt,1 vt(i1) = t_min+dt*dfloat(i1) 100 continue C de = (e_max-e_min)/dfloat(ne) do 110 i2=0,ne,1 ve(i2) = e_min+de*dfloat(i2) 110 continue C do 120 i1=0,nt,1 do 130 i2=0,ne,1 read(file_alm,*) fac(i1,i2) do 140 l=0,tl_max,1 do 150 m=(-l),l,1 read(file_alm,*) alm(i1,i2,l,m) 150 continue 140 continue 130 continue 120 continue C n = dim nc = 0 nd = 0 x(1) = (e_min+e_max)/2.d00 x(2) = (t_min+t_max)/2.d00 x(3) = 0.d00 x(4) = 65.4321d00 x(5) = 123.456d00 c d(1) = de c d(2) = dt d(1) = (e_max-e_min)/5.d00 d(2) = (t_max-t_min)/5.d00 d(3) = 180.d00 d(4) = 90.d00 d(5) = 180.d00 wi = 1.d-30 C do 200 k1=1,nev,1 do 210 k2=1,20,1 C ++++++++++++++++++++++++++++++++++++ call metropolis(n,nc,nd,wi,d,x,work) C ++++++++++++++++++++++++++++++++++++ if (x(3).lt.0.d00) x(3)=x(3)+360.d00 if (x(3).gt.360.d00) x(3)=x(3)-360.d00 if (x(5).lt.0.d00) x(5)=x(5)+360.d00 if (x(5).gt.360.d00) x(5)=x(5)-360.d00 210 continue c write(101,2000) x(1),x(2),x(3),x(4),x(5) c write(102,*) k1,dfloat(nd)/dfloat(nc) e_rel = x(1) theta_cm = x(2) phi_cm = x(3) theta_rel = x(4) phi_rel = x(5) C +++++++++++++++++++++++++++++++++++++++++++++++++ call rkin(e_rel,theta_cm,phi_cm,theta_rel,phi_rel & ,pb,pc,px) C +++++++++++++++++++++++++++++++++++++++++++++++++ write(file_ev,2001) pb,pc,px c write(103,*) (pb(2)+pc(2)+px(2)) c & ,(pb(3)+pc(3)+px(3)) c & ,(pb(4)+pc(4)+px(4)) 200 continue C acc = dfloat(nd)/dfloat(nc) C 2000 format(5(2x,f12.8)) 2001 format(12(1x,e14.7)) C return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ******************************************************* subroutine rkin(e_rel,theta_cm,phi_cm,theta_rel,phi_rel & ,pb,pc,px) C ******************************************************* implicit none integer*4 i1,i2 real*8 e_rel,theta_cm,phi_cm,theta_rel,phi_rel & ,pb(4),pc(4),px(4),tmp1(4),tmp2(4),tmp3(4) & ,pi,tpi,radfac,m_x,m_a,m_b,m_c,e_i,t_p & ,m_aa,s1,s2,p2,p3,pa,eaa,vaa,v1(3),v2(3) & ,thcm,phcm,thr,phr,g1,g2 & ,n1(3),n2(3),n3(3),lt1(4,4),lt2(4,4) C common /simpar0/ pi,tpi,radfac common /simpar1/ m_x,m_a,m_b,m_c,e_i,t_p C s1 = (m_a+m_x)*(m_a+m_x)+2.d00*t_p*m_x m_aa = m_b+m_c+e_rel s2 = m_aa*m_aa p2 = s1-s2-m_x*m_x p2 = dsqrt((p2*p2-4.d00*s2*m_x*m_x)/(4.d00*s1)) p3 = s2-m_b*m_b-m_c*m_c p3 = dsqrt((p3*p3-4.d00*m_b*m_b*m_c*m_c)/(4.d00*s2)) C pa = dsqrt(t_p*(2.d00*m_a+t_p)) n1(1) = 0.d00 n1(2) = 0.d00 n1(3) = 1.d00 C velocity of a+x cm system in lab system v1(1) = 0.d00 v1(2) = 0.d00 v1(3) = pa/(m_a+m_x+t_p) g1 = 1.d00/dsqrt(1.d00-v1(3)*v1(3)) C C Lorentz transformation from a+x cm system to lab system lt1(1,1) = g1 lt1(1,2) = g1*v1(1) lt1(1,3) = g1*v1(2) lt1(1,4) = g1*v1(3) lt1(2,1) = lt1(1,2) lt1(2,2) = 1.d00+(g1-1.d00)*n1(1)*n1(1) lt1(2,3) = (g1-1.d00)*n1(1)*n1(2) lt1(2,4) = (g1-1.d00)*n1(1)*n1(3) lt1(3,1) = lt1(1,3) lt1(3,2) = lt1(2,3) lt1(3,3) = 1.d00+(g1-1.d00)*n1(2)*n1(2) lt1(3,4) = (g1-1.d00)*n1(2)*n1(3) lt1(4,1) = lt1(1,4) lt1(4,2) = lt1(2,4) lt1(4,3) = lt1(3,4) lt1(4,4) = 1.d00+(g1-1.d00)*n1(3)*n1(3) C thcm = theta_cm*radfac phcm = phi_cm*radfac C n2(1) = dsin(thcm)*dcos(phcm) n2(2) = dsin(thcm)*dsin(phcm) n2(3) = dcos(thcm) C eaa = dsqrt(m_aa*m_aa+p2*p2) C velocity of A=b+c cm system in a+x cm system vaa = p2/eaa v2(1) = vaa*n2(1) v2(2) = vaa*n2(2) v2(3) = vaa*n2(3) g2 = 1.d00/dsqrt(1.d00-vaa*vaa) C C Lorentz transformation from A=b+c cm system to a+x cm system lt2(1,1) = g2 lt2(1,2) = g2*v2(1) lt2(1,3) = g2*v2(2) lt2(1,4) = g2*v2(3) lt2(2,1) = lt2(1,2) lt2(2,2) = 1.d00+(g2-1.d00)*n2(1)*n2(1) lt2(2,3) = (g2-1.d00)*n2(1)*n2(2) lt2(2,4) = (g2-1.d00)*n2(1)*n2(3) lt2(3,1) = lt2(1,3) lt2(3,2) = lt2(2,3) lt2(3,3) = 1.d00+(g2-1.d00)*n2(2)*n2(2) lt2(3,4) = (g2-1.d00)*n2(2)*n2(3) lt2(4,1) = lt2(1,4) lt2(4,2) = lt2(2,4) lt2(4,3) = lt2(3,4) lt2(4,4) = 1.d00+(g2-1.d00)*n2(3)*n2(3) c do 200 i1=1,4,1 c write(200,1000) (lt2(i1,i2),i2=1,4) c 200 continue c stop thr = theta_rel*radfac phr = (phi_cm+phi_rel)*radfac C n3(1) = dsin(thr)*dcos(phr) n3(2) = dsin(thr)*dsin(phr) n3(3) = dcos(thr) C C four momentum of nucleus b in A=b+c cm system tmp1(1) = dsqrt(m_b*m_b+p3*p3) tmp1(2) = p3*n3(1) tmp1(3) = p3*n3(2) tmp1(4) = p3*n3(3) C C four momentum of nucleus c in A=b+c cm system tmp2(1) = dsqrt(m_c*m_c+p3*p3) tmp2(2) = -tmp1(2) tmp2(3) = -tmp1(3) tmp2(4) = -tmp1(4) C C four momentum of nucleus b in a+x cm system do 100 i1=1,4,1 tmp3(i1) = 0.d00 do 110 i2=1,4,1 tmp3(i1) = tmp3(i1)+lt2(i1,i2)*tmp1(i2) 110 continue 100 continue C four momentum of nucleus b in lab system do 120 i1=1,4,1 pb(i1) = 0.d00 do 130 i2=1,4,1 pb(i1) = pb(i1)+lt1(i1,i2)*tmp3(i2) 130 continue 120 continue C C four momentum of nucleus c in a+x cm system do 140 i1=1,4,1 tmp3(i1) = 0.d00 c tmp3(i1) = tmp2(i1) do 150 i2=1,4,1 tmp3(i1) = tmp3(i1)+lt2(i1,i2)*tmp2(i2) 150 continue 140 continue C four momentum of nucleus c in lab system do 160 i1=1,4,1 pc(i1) = 0.d00 c pc(i1) = tmp3(i1) do 170 i2=1,4,1 pc(i1) = pc(i1)+lt1(i1,i2)*tmp3(i2) 170 continue 160 continue C C four momentum of nucleus x in lab system px(2) = -pb(2)-pc(2) px(3) = -pb(3)-pc(3) px(4) = pa-pb(4)-pc(4) px(1) = dsqrt(m_x*m_x+px(2)*px(2)+px(3)*px(3)+px(4)*px(4)) C return end C ******************************************* subroutine metropolis(n,nc,nd,wi,d,xi,work) C ******************************************* C written by Stefan Typel C version 1.1 C 2001/04/15 C C wi must be set to unequal 0.d00 before the first call of metropolis C nc and nd must be set to 0 before the first call of metropolis C no correction for correlation C implicit none integer*4 n,k,nc,nd real*4 rand real*8 prob,r,wi,wo,xi(5),work(5),d(5),rr C nc=nc+1 do 100 k=1,n,1 work(k)=xi(k)+d(k)*(2.d00*rand(0)-1.d00) 100 continue wo=prob(n,work) 300 continue rr = rand(0) if (rr.eq.0.d00) goto 300 r=wo/wi if (rr.gt.r) goto 200 nd=nd+1 wi=wo do 110 k=1,n,1 xi(k)=work(k) 110 continue 200 continue C return end C *********************************** double precision function prob(n,x) C *********************************** implicit none integer*4 n,nt,ne real*8 x(5),e_rel,theta_cm,theta_rel,phi_rel,d3s & ,t_min,t_max,dt,e_min,e_max,de & ,pi,tpi,radfac C common /simpar0/ pi,tpi,radfac common /simpar2/ t_min,t_max,dt,e_min,e_max,de,nt,ne C prob=0.d00 C e_rel = x(1) if (e_rel.gt.e_max) return if (e_rel.lt.e_min) return theta_cm = x(2) if (theta_cm.gt.t_max) return if (theta_cm.lt.t_min) return theta_rel = x(4) if (theta_rel.gt.180.d00) return if (theta_rel.lt.0.d00) return phi_rel = x(5) C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ call get_d3sigma(theta_cm,e_rel,theta_rel,phi_rel,d3s) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++ C prob = d3s & *dabs(dsin(theta_cm*radfac)*dsin(theta_rel*radfac)) C return end C ************************************************************ subroutine get_d3sigma(theta_cm,e_rel,theta_rel,phi_rel,d3s) C ************************************************************ implicit none integer*4 nt,ne,nt_max,ne_max,l_max,tl_max,i1,i2,l,m parameter (nt_max=90,ne_max=501,l_max=3,tl_max=2*l_max) real*8 theta_cm,e_rel,theta_rel,phi_rel,pe,pt & ,t_min,t_max,dt,e_min,e_max,de & ,vt(0:nt_max),ve(0:ne_max),fac(0:nt_max,0:ne_max) complex*16 ylm(0:tl_max,-tl_max:tl_max) & ,d3s,d3s00,d3s01,d3s10,d3s11,cd3s & ,alm(0:nt_max,0:ne_max,0:tl_max,-tl_max:tl_max) C common /simpar2/ t_min,t_max,dt,e_min,e_max,de,nt,ne common /simpar3/ vt,ve common /simpar4/ alm common /simpar5/ fac common /ylmtp/ ylm C i1 = idnint((theta_cm-t_min)/dt-0.5d00) if (i1.lt.0) i1=0 if (i1.ge.nt) i1=nt-1 pt = (theta_cm-vt(i1))/dt C i2 = idnint((e_rel-e_min)/de-0.5d00) if (i2.lt.0) i2=0 if (i2.ge.ne) i2=ne-1 pe = (e_rel-ve(i2))/de C C +++++++++++++++++++++++++++++++ call get_ylm(theta_rel,phi_rel) C +++++++++++++++++++++++++++++++ C d3s00=0.d00 d3s10=0.d00 d3s01=0.d00 d3s11=0.d00 C do 100 l=0,tl_max,1 do 110 m=(-l),l,1 d3s00=d3s00+alm(i1 ,i2 ,l,m)*ylm(l,m) d3s10=d3s10+alm(i1+1,i2 ,l,m)*ylm(l,m) d3s01=d3s01+alm(i1 ,i2+1,l,m)*ylm(l,m) d3s11=d3s11+alm(i1+1,i2+1,l,m)*ylm(l,m) 110 continue 100 continue C d3s00 = d3s00*fac(i1 ,i2 ) d3s10 = d3s10*fac(i1+1,i2 ) d3s01 = d3s01*fac(i1 ,i2+1) d3s11 = d3s11*fac(i1+1,i2+1) C cd3s = (1.d00-pt)*((1.d00-pe)*d3s00+pe*d3s01) & +pt*((1.d00-pe)*d3s10+pe*d3s11) C d3s = dabs(dreal(cd3s)) C return end