created branch for code restructuring (progressive removal of common statements, dynamic allocation, etc). curr_int subroutine (and dependencies) partially updated.
This commit is contained in:
parent
1df304f221
commit
3b24d5d58d
19
Makefile
19
Makefile
@ -2,7 +2,8 @@
|
||||
EXE=gray
|
||||
|
||||
# Objects list
|
||||
OBJ=gray.o grayl.o reflections.o green_func_p.o \
|
||||
MAINOBJ=gray.o
|
||||
OTHOBJ=dqagmv.o grayl.o reflections.o green_func_p.o \
|
||||
const_and_precisions.o itm_constants.o itm_types.o
|
||||
|
||||
# Alternative search paths
|
||||
@ -11,32 +12,32 @@ vpath %.f src
|
||||
|
||||
# Fortran compiler name and flags
|
||||
FC=gfortran
|
||||
FFLAGS=-Wall -fcheck=all
|
||||
FFLAGS=-O3 #-Wall -g -fcheck=all
|
||||
|
||||
DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
|
||||
|
||||
all: $(EXE)
|
||||
|
||||
# Build executable from object files
|
||||
$(EXE): $(OBJ)
|
||||
$(EXE): $(MAINOBJ) $(OTHOBJ)
|
||||
$(FC) $(FFLAGS) -o $@ $^
|
||||
|
||||
# Dependencies on modules
|
||||
gray.o: green_func_p.o reflections.o
|
||||
gray.o: dqagmv.o green_func_p.o reflections.o
|
||||
green_func_p.o: const_and_precisions.o
|
||||
const_and_precisions.o: itm_types.o itm_constants.o
|
||||
itm_constants.o: itm_types.o
|
||||
|
||||
# General object compilation command
|
||||
%.o: %.f90
|
||||
$(FC) $(FFLAGS) -c $<
|
||||
$(FC) $(FFLAGS) -c $^
|
||||
|
||||
gray.o:gray.f green_func_p.o
|
||||
%.o: %.f
|
||||
$(FC) $(FFLAGS) -c $^
|
||||
|
||||
gray.o:gray.f
|
||||
$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $<
|
||||
|
||||
grayl.o:grayl.f
|
||||
$(FC) $(FFLAGS) -c $^
|
||||
|
||||
.PHONY: clean install
|
||||
# Remove output files
|
||||
clean:
|
||||
|
||||
1695
src/dqagmv.f
Normal file
1695
src/dqagmv.f
Normal file
File diff suppressed because it is too large
Load Diff
277
src/gray.f
277
src/gray.f
@ -5798,10 +5798,10 @@ c
|
||||
implicit none
|
||||
real*8 anum,denom,resp,resj,effjcd,coullog,dens,tekev
|
||||
real*8 anucc,canucc,ddens
|
||||
real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff
|
||||
real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff,psinv
|
||||
real*8 rbn,rbx,btot,bmin,bmax,alams,fp0s,pa,fc
|
||||
real*8 fjch,fjncl,fjch0,fconic
|
||||
real*8 cst,cst2
|
||||
real*8 cst,cst2,eccdpar(5)
|
||||
integer ieccd,nharm,nhf,nhn
|
||||
external fjch,fjncl,fjch0
|
||||
|
||||
@ -5818,11 +5818,10 @@ c
|
||||
common/tete/tekev
|
||||
common/dens/dens,ddens
|
||||
common/zz/Zeff
|
||||
common/psival/psinv
|
||||
common/btot/btot
|
||||
common/bmxmn/bmax,bmin
|
||||
common/fc/fc
|
||||
common/ncl/rbx
|
||||
common/cohen/rbn,alams,pa,fp0s
|
||||
common/cst/cst,cst2
|
||||
|
||||
anum=0.0d0
|
||||
@ -5833,6 +5832,8 @@ c
|
||||
anucc=canucc*dens*coullog
|
||||
|
||||
c nhf=nharm
|
||||
|
||||
eccdpar(1)=zeff
|
||||
|
||||
select case(ieccd)
|
||||
|
||||
@ -5852,8 +5853,12 @@ c fp0s= P_a (alams)
|
||||
alams=sqrt(1.0d0-bmin/bmax)
|
||||
pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0
|
||||
fp0s=fconic(alams,pa,0)
|
||||
eccdpar(2)=rbn
|
||||
eccdpar(3)=alams
|
||||
eccdpar(4)=pa
|
||||
eccdpar(5)=fp0s
|
||||
do nhn=nharm,nhf
|
||||
call curr_int(fjch,resj,resp)
|
||||
call curr_int(fjch,eccdpar,5,resj,resp)
|
||||
anum=anum+resj
|
||||
denom=denom+resp
|
||||
end do
|
||||
@ -5862,7 +5867,7 @@ c fp0s= P_a (alams)
|
||||
cst=0.0d0
|
||||
cst2=0.0d0
|
||||
do nhn=nharm,nhf
|
||||
call curr_int(fjch0,resj,resp)
|
||||
call curr_int(fjch0,eccdpar,1,resj,resp)
|
||||
anum=anum+resj
|
||||
denom=denom+resp
|
||||
end do
|
||||
@ -5877,8 +5882,12 @@ c rzfc=(1+Zeff)/fc
|
||||
if(cst2.lt.1d-6) cst2=0.0d0
|
||||
cst=sqrt(cst2)
|
||||
call SpitzFuncCoeff(Tekev,Zeff,fc)
|
||||
eccdpar(2) = tekev
|
||||
eccdpar(3) = fc
|
||||
eccdpar(4) = rbx
|
||||
eccdpar(5) = psinv
|
||||
do nhn=nharm,nhf
|
||||
call curr_int(fjncl,resj,resp)
|
||||
call curr_int(fjncl,eccdpar,5,resj,resp)
|
||||
anum=anum+resj
|
||||
denom=denom+resp
|
||||
end do
|
||||
@ -5897,23 +5906,30 @@ c
|
||||
c
|
||||
c
|
||||
c
|
||||
subroutine curr_int(fcur,resj,resp)
|
||||
subroutine curr_int(fcur,eccdpar,necp,resj,resp)
|
||||
implicit real*8(a-h,o-z)
|
||||
parameter(epsa=0.0d0,epsr=1.0d-2)
|
||||
parameter(xxcr=16.0d0)
|
||||
parameter (lw=5000,liw=lw/4)
|
||||
dimension w(lw),iw(liw)
|
||||
parameter (npar=20)
|
||||
complex*16 ex,ey,ez
|
||||
dimension w(lw),iw(liw),eccdpar(necp),apar(npar)
|
||||
external fcur,fpp
|
||||
|
||||
common/nhn/nhn
|
||||
common/ygyg/yg
|
||||
common/nplr/anpl,anpr
|
||||
common/amut/amu
|
||||
common/gg/uplp,uplm,ygn
|
||||
common/ierr/ierr
|
||||
common/iokh/iokhawa
|
||||
common/cst/cst,cst2
|
||||
|
||||
c used and passed
|
||||
common/ygyg/yg
|
||||
common/nplr/anpl,anpr
|
||||
common/amut/amu
|
||||
common/nhn/nhn
|
||||
c passed unused
|
||||
common/epolar/ex,ey,ez
|
||||
common/nprw/anprre,anprim
|
||||
common/ithn/ithn
|
||||
|
||||
c EC power and current densities
|
||||
|
||||
anpl2=anpl*anpl
|
||||
@ -5945,10 +5961,30 @@ c
|
||||
if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl
|
||||
if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
|
||||
duu=abs(uu1-uu2)
|
||||
c
|
||||
apar(1) = yg
|
||||
apar(2) = anpl
|
||||
apar(3) = amu
|
||||
apar(4) = anprre
|
||||
apar(5) = dble(ex)
|
||||
apar(6) = dimag(ex)
|
||||
apar(7) = dble(ey)
|
||||
apar(8) = dimag(ey)
|
||||
apar(9) = dble(ez)
|
||||
apar(10) = dimag(ez)
|
||||
apar(11) = dble(ithn)
|
||||
apar(12) = dble(nhn)
|
||||
apar(13) = uplp
|
||||
apar(14) = uplm
|
||||
apar(15) = ygn
|
||||
c
|
||||
do i=1,necp
|
||||
apar(15+i) = eccdpar(i)
|
||||
end do
|
||||
c
|
||||
if(duu.gt.1.d-6) then
|
||||
call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ier,
|
||||
. liw,lw,last,iw,w)
|
||||
call dqagsmv(fpp,uu1,uu2,apar,npar,epsa,epsr,resp,epp,neval,
|
||||
. ier,liw,lw,last,iw,w)
|
||||
if (ier.gt.0) ierr=90
|
||||
end if
|
||||
c
|
||||
@ -5960,7 +5996,7 @@ c resonance curve does not cross the trapping region
|
||||
c
|
||||
iokhawa=0
|
||||
if(duu.gt.1.d-4) then
|
||||
call dqags(fcur,uu1,uu2,epsa,epsr,
|
||||
call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
|
||||
. resj,ej,neval,ier,liw,lw,last,iw,w)
|
||||
if (ier.gt.0) ierr=91
|
||||
end if
|
||||
@ -5982,7 +6018,7 @@ c
|
||||
if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
|
||||
duu=abs(uu1-uu2)
|
||||
if(duu.gt.1.d-6) then
|
||||
call dqags(fcur,uu1,uu2,epsa,epsr,
|
||||
call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
|
||||
. resj1,ej1,neval,ier,liw,lw,last,iw,w)
|
||||
if (ier.gt.0) then
|
||||
if (abs(resj1).lt.1.0d-10) then
|
||||
@ -6003,7 +6039,7 @@ c
|
||||
if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
|
||||
duu=abs(uu1-uu2)
|
||||
if(duu.gt.1.d-6) then
|
||||
call dqags(fcur,uu1,uu2,epsa,epsr,
|
||||
call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
|
||||
. resj2,ej2,neval,ier,liw,lw,last,iw,w)
|
||||
if (ier.gt.0) then
|
||||
if(ier.ne.2) ierr=93
|
||||
@ -6023,20 +6059,45 @@ c ith=0 : polarization term = const
|
||||
c ith=1 : polarization term Larmor radius expansion to lowest order
|
||||
c ith=2 : full polarization term (J Bessel)
|
||||
c
|
||||
function fpp(upl)
|
||||
function fpp(upl,extrapar,npar)
|
||||
c
|
||||
c integration variable upl passed explicitly. other variables passed
|
||||
c as array of extra parameters of length npar=size(extrapar)
|
||||
c
|
||||
c extrapar(1) = yg
|
||||
c extrapar(2) = anpl
|
||||
c extrapar(3) = amu
|
||||
c extrapar(4) = Re(anprw)
|
||||
c extrapar(5) = Re(ex)
|
||||
c extrapar(6) = Im(ex)
|
||||
c extrapar(7) = Re(ey)
|
||||
c extrapar(8) = Im(ey)
|
||||
c extrapar(9) = Re(ez)
|
||||
c extrapar(10) = Im(ez)
|
||||
c extrapar(11) = double(ithn)
|
||||
c extrapar(12) = double(nhn)
|
||||
c extrapar(13) = uplp
|
||||
c extrapar(14) = uplm
|
||||
c extrapar(15) = ygn
|
||||
c
|
||||
implicit real*8 (a-h,o-z)
|
||||
complex*16 ex,ey,ez,ui,emxy,epxy
|
||||
parameter(ui=(0.0d0,1.0d0))
|
||||
|
||||
common/ygyg/yg
|
||||
common/nplr/anpl,anpr
|
||||
common/amut/amu
|
||||
common/gg/uplp,uplm,ygn
|
||||
common/epolar/ex,ey,ez
|
||||
common/nprw/anprre,anprim
|
||||
common/ithn/ithn
|
||||
common/nhn/nhn
|
||||
dimension extrapar(npar)
|
||||
|
||||
yg=extrapar(1)
|
||||
anpl=extrapar(2)
|
||||
amu=extrapar(3)
|
||||
anprre=extrapar(4)
|
||||
ex=dcmplx(extrapar(5),extrapar(6))
|
||||
ey=dcmplx(extrapar(7),extrapar(8))
|
||||
ez=dcmplx(extrapar(9),extrapar(10))
|
||||
ithn=int(extrapar(11))
|
||||
nhn=int(extrapar(12))
|
||||
uplp=extrapar(13)
|
||||
uplm=extrapar(14)
|
||||
ygn=extrapar(15)
|
||||
|
||||
upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm)
|
||||
gam=anpl*upl+ygn
|
||||
ee=exp(-amu*(gam-1))
|
||||
@ -6076,15 +6137,57 @@ c fjch integrand for Cohen model with trapping
|
||||
c fjch0 integrand for Cohen model without trapping
|
||||
c fjncl integrand for momentum conserv. model K(u) from Maruschenko
|
||||
c gg=F(u)/u with F(u) as in Cohen paper
|
||||
|
||||
function fjch(upl)
|
||||
c
|
||||
function fjch(upl,extrapar,npar)
|
||||
c
|
||||
c integration variable upl passed explicitly. Other variables passed
|
||||
c as array of extra parameters of length npar=size(extrapar).
|
||||
c variables with index 1..15 must be passed to fpp
|
||||
c variable with index 16 is zeff
|
||||
c variables with index gt 16 are specific of the cd model
|
||||
c
|
||||
c extrapar(1) = yg
|
||||
c extrapar(2) = anpl
|
||||
c extrapar(3) = amu
|
||||
c extrapar(4) = Re(anprw)
|
||||
c extrapar(5) = Re(ex)
|
||||
c extrapar(6) = Im(ex)
|
||||
c extrapar(7) = Re(ey)
|
||||
c extrapar(8) = Im(ey)
|
||||
c extrapar(9) = Re(ez)
|
||||
c extrapar(10) = Im(ez)
|
||||
c extrapar(11) = double(ithn)
|
||||
c extrapar(12) = double(nhn)
|
||||
c extrapar(13) = uplp
|
||||
c extrapar(14) = uplm
|
||||
c extrapar(15) = ygn
|
||||
c
|
||||
c extrapar(16) = zeff
|
||||
c extrapar(17) = rb
|
||||
c extrapar(18) = alams
|
||||
c extrapar(19) = pa
|
||||
c extrapar(20) = fp0s
|
||||
c
|
||||
implicit real*8 (a-h,o-z)
|
||||
complex*16 ex,ey,ez
|
||||
dimension extrapar(npar)
|
||||
|
||||
anpl=extrapar(2)
|
||||
anprre=extrapar(4)
|
||||
ex=dcmplx(extrapar(5),extrapar(6))
|
||||
ey=dcmplx(extrapar(7),extrapar(8))
|
||||
ez=dcmplx(extrapar(9),extrapar(10))
|
||||
ithn=int(extrapar(11))
|
||||
nhn=int(extrapar(12))
|
||||
uplp=extrapar(13)
|
||||
uplm=extrapar(14)
|
||||
ygn=extrapar(15)
|
||||
zeff=extrapar(16)
|
||||
rb=extrapar(17)
|
||||
alams=extrapar(18)
|
||||
pa=extrapar(19)
|
||||
fp0s=extrapar(20)
|
||||
|
||||
common/gg/uplp,uplm,ygn
|
||||
common/nplr/anpl,anpr
|
||||
common/zz/Zeff
|
||||
common/cohen/rb,alams,pa,fp0s
|
||||
|
||||
upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm)
|
||||
gam=anpl*upl+ygn
|
||||
u2=gam*gam-1.0d0
|
||||
@ -6111,17 +6214,52 @@ c gg=F(u)/u with F(u) as in Cohen paper
|
||||
eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam)
|
||||
|
||||
if(upl.lt.0.0d0) eta=-eta
|
||||
fjch=eta*fpp(upl)
|
||||
fjch=eta*fpp(upl,extrapar,npar)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
function fjch0(upl)
|
||||
function fjch0(upl,extrapar,npar)
|
||||
c
|
||||
c integration variable upl passed explicitly. Other variables passed
|
||||
c as array of extra parameters of length npar=size(extrapar).
|
||||
c variables with index 1..15 must be passed to fpp
|
||||
c variable with index 16 is zeff
|
||||
c variables with index gt 16 are specific of the cd model
|
||||
c
|
||||
c extrapar(1) = yg
|
||||
c extrapar(2) = anpl
|
||||
c extrapar(3) = amu
|
||||
c extrapar(4) = Re(anprw)
|
||||
c extrapar(5) = Re(ex)
|
||||
c extrapar(6) = Im(ex)
|
||||
c extrapar(7) = Re(ey)
|
||||
c extrapar(8) = Im(ey)
|
||||
c extrapar(9) = Re(ez)
|
||||
c extrapar(10) = Im(ez)
|
||||
c extrapar(11) = double(ithn)
|
||||
c extrapar(12) = double(nhn)
|
||||
c extrapar(13) = uplp
|
||||
c extrapar(14) = uplm
|
||||
c extrapar(15) = ygn
|
||||
c
|
||||
c extrapar(16) = zeff
|
||||
c
|
||||
implicit real*8 (a-h,o-z)
|
||||
common/gg/uplp,uplm,ygn
|
||||
common/nplr/anpl,anpr
|
||||
common/zz/Zeff
|
||||
complex*16 ex,ey,ez
|
||||
dimension extrapar(npar)
|
||||
|
||||
anpl=extrapar(2)
|
||||
anprre=extrapar(4)
|
||||
ex=dcmplx(extrapar(5),extrapar(6))
|
||||
ey=dcmplx(extrapar(7),extrapar(8))
|
||||
ez=dcmplx(extrapar(9),extrapar(10))
|
||||
ithn=int(extrapar(11))
|
||||
nhn=int(extrapar(12))
|
||||
ygn=extrapar(15)
|
||||
zeff=extrapar(16)
|
||||
|
||||
gam=anpl*upl+ygn
|
||||
u2=gam*gam-1.0d0
|
||||
u=sqrt(u2)
|
||||
@ -6135,24 +6273,54 @@ c gg=F(u)/u with F(u) as in Cohen paper
|
||||
gg=u*gu/z5
|
||||
dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5
|
||||
eta=anpl*gg+gam*upl*dgg/u
|
||||
fjch0=eta*fpp(upl)
|
||||
fjch0=eta*fpp(upl,extrapar,npar)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
function fjncl(upl)
|
||||
use green_func_p
|
||||
function fjncl(upl,extrapar,npar)
|
||||
c
|
||||
c integration variable upl passed explicitly. Other variables passed
|
||||
c as array of extra parameters of length npar=size(extrapar).
|
||||
c variables with index 1..15 must be passed to fpp
|
||||
c variable with index 16 is zeff
|
||||
c variables with index gt 16 are specific of the cd model
|
||||
c
|
||||
c extrapar(1) = yg
|
||||
c extrapar(2) = anpl
|
||||
c extrapar(3) = amu
|
||||
c extrapar(4) = Re(anprw)
|
||||
c extrapar(5) = Re(ex)
|
||||
c extrapar(6) = Im(ex)
|
||||
c extrapar(7) = Re(ey)
|
||||
c extrapar(8) = Im(ey)
|
||||
c extrapar(9) = Re(ez)
|
||||
c extrapar(10) = Im(ez)
|
||||
c extrapar(11) = double(ithn)
|
||||
c extrapar(12) = double(nhn)
|
||||
c extrapar(13) = uplp
|
||||
c extrapar(14) = uplm
|
||||
c extrapar(15) = ygn
|
||||
c
|
||||
c extrapar(16) = zeff
|
||||
c extrapar(17) = tekev
|
||||
c extrapar(18) = fc
|
||||
c extrapar(19) = hb
|
||||
c extrapar(20) = psin
|
||||
c
|
||||
use green_func_p, only: GenSpitzFunc
|
||||
implicit real*8 (a-h,o-z)
|
||||
dimension extrapar(npar)
|
||||
|
||||
common/gg/uplp,uplm,ygn
|
||||
common/nplr/anpl,anpr
|
||||
common/fc/fc
|
||||
common/ncl/hb
|
||||
common/psival/psinv
|
||||
common/amut/amu
|
||||
common/tete/tekev
|
||||
common/zz/Zeff
|
||||
anpl=extrapar(2)
|
||||
amu=extrapar(3)
|
||||
ygn=extrapar(15)
|
||||
zeff=extrapar(16)
|
||||
tekev=extrapar(17)
|
||||
fc=extrapar(18)
|
||||
hb=extrapar(19)
|
||||
psin=extrapar(20)
|
||||
|
||||
gam=anpl*upl+ygn
|
||||
u2=gam*gam-1.0d0
|
||||
@ -6165,12 +6333,11 @@ c gg=F(u)/u with F(u) as in Cohen paper
|
||||
dfk=dfk*(2.0d0/amu)*bth
|
||||
|
||||
alam=upr2/u2/hb
|
||||
psi=psinv
|
||||
call vlambda(alam,psi,fu,dfu)
|
||||
call vlambda(alam,psin,fu,dfu)
|
||||
|
||||
eta=gam*fu*dfk/u-2.0d0*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb
|
||||
if(upl.lt.0) eta=-eta
|
||||
fjncl=eta*fpp(upl)
|
||||
fjncl=eta*fpp(upl,extrapar,npar)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user