istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src: remove unnecessary one, zero uses

Browse Source
This commit is contained in:
Michele Guerini Rocco 2024-09-23 22:16:33 +02:00 committed by rnhmjoj
parent 80782a58fc
commit 24e0e6e472
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
7 changed files with 221 additions and 250 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
src/dispersion.f90
View File

@ -145,7 +145,7 @@ pure subroutine harmnumber(Y, mu, Npl2, weakly, nhmin, nhmax)
! Yc = 1 -½ N∥² (weakly relativistic) ! Yc = 1 -½ N∥² (weakly relativistic)
! !
if (weakly) then if (weakly) then
Yc = max(one - npl2/2, zero) Yc = max(1 - npl2/2, zero)
else else
Yc = sqrt(max(1 - npl2, zero)) Yc = sqrt(max(1 - npl2, zero))
end if end if
@ -196,7 +196,7 @@ pure subroutine harmnumber(Y, mu, Npl2, weakly, nhmin, nhmax)
else else
rdu2 = Yn**2 - Yc**2 rdu2 = Yn**2 - Yc**2
gamma = (Yn - sqrt(Npl2*rdu2))/(1 - Npl2) gamma = (Yn - sqrt(Npl2*rdu2))/(1 - Npl2)
argexp = mu*(gamma - one) argexp = mu*(gamma - 1)
end if end if
if (argexp <= expcr) then if (argexp <= expcr) then
@ -406,7 +406,7 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, &
Npr2 = Npr2_prev + 0.05_wp_ * (Npr2 - Npr2_prev) / abs(Npr2 - Npr2_prev) Npr2 = Npr2_prev + 0.05_wp_ * (Npr2 - Npr2_prev) / abs(Npr2 - Npr2_prev)
! Again, make sure that we have a damped EM wave and not a ! Again, make sure that we have a damped EM wave and not a
! Bernstein-like wave (see above) ! Bernstein-like wave (see above)
if (real(sqrt(Npr2)) * aimag(sqrt(Npr2)) < zero) & if (real(sqrt(Npr2)) * aimag(sqrt(Npr2)) < 0) &
Npr2 = conjg(sqrt(Npr2))**2 Npr2 = conjg(sqrt(Npr2))**2
end if end if
end block modify_fixed_point end block modify_fixed_point
@ -464,7 +464,7 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, &
if (sox < 0) then if (sox < 0) then
e(3) = 1 e(3) = 1
else else
e(1) = sqrt(one/(1 + abs(-eps(1,1)/eps(1,2))**2)) e(1) = sqrt(1/(1 + abs(-eps(1,1)/eps(1,2))**2))
e(2) = -e(1)*eps(1,1)/eps(1,2) e(2) = -e(1)*eps(1,1)/eps(1,2)
end if end if
end if end if
@ -673,7 +673,7 @@ subroutine hermitian(rr,yg,mu,npl,cr,fast,lrm)
do n=-lrm,lrm do n=-lrm,lrm
do k=0,2 do k=0,2
do m=0,lrm do m=0,lrm
rr(n,k,m)=zero rr(n,k,m) = 0
end do end do
end do end do
end do end do
@ -893,7 +893,7 @@ subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm,error)
do n=-lrm,lrm do n=-lrm,lrm
do k=0,2 do k=0,2
do m=0,lrm do m=0,lrm
rr(n,k,m)=zero rr(n,k,m) = 0
end do end do
end do end do
end do end do
@ -1059,7 +1059,7 @@ function fhermit(t,apar,npar)
zm2=zm*zm zm2=zm*zm
zm3=zm2*zm zm3=zm2*zm
call calcei3(zm,fe0m) call calcei3(zm,fe0m)
ffe=zero ffe = 0
uplh=upl**ih uplh=upl**ih
if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2
if(m.eq.1) ffe=(one+s*(one-zm*fe0m))*uplh/mu2 if(m.eq.1) ffe=(one+s*(one-zm*fe0m))*uplh/mu2
@ -1089,7 +1089,7 @@ subroutine antihermitian(ri,yg,mu,npl,ci,lrm)
do n=1,lrm do n=1,lrm
do k=0,2 do k=0,2
do m=1,lrm do m=1,lrm
ri(n,k,m)=zero ri(n,k,m) = 0
end do end do
end do end do
end do end do
@ -1101,7 +1101,7 @@ subroutine antihermitian(ri,yg,mu,npl,ci,lrm)
do n=1,lrm do n=1,lrm
ygn=n*yg ygn=n*yg
rdu2=ygn**2-dnl rdu2=ygn**2-dnl
if(rdu2.gt.zero) then if(rdu2.gt.0) then
rdu=sqrt(rdu2) rdu=sqrt(rdu2)
du=rdu/dnl du=rdu/dnl
ub=npl*ygn/dnl ub=npl*ygn/dnl
@ -1231,17 +1231,17 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error)
phim=sqrt(abs(phi2)) phim=sqrt(abs(phi2))
if (alpha.ge.0) then if (alpha.ge.0) then
xp=psi-phim xp=psi-phim
yp=zero yp=0
xm=-psi-phim xm=-psi-phim
ym=zero ym=0
x0=-phim x0=-phim
y0=zero y0=0
else else
xp=psi xp=psi
yp=phim yp=phim
xm=-psi xm=-psi
ym=phim ym=phim
x0=zero x0=0
y0=phim y0=phim
end if end if
call zetac (xp,yp,zrp,zip,iflag) call zetac (xp,yp,zrp,zip,iflag)
@ -1295,19 +1295,19 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error)
alpha=npl*npl/2.0_wp_-is*yg-one alpha=npl*npl/2.0_wp_-is*yg-one
phi2=mu*alpha phi2=mu*alpha
phim=sqrt(abs(phi2)) phim=sqrt(abs(phi2))
if (alpha.ge.zero) then if (alpha.ge.0) then
xp=psi-phim xp=psi-phim
yp=zero yp=0
xm=-psi-phim xm=-psi-phim
ym=zero ym=0
x0=-phim x0=-phim
y0=zero y0=0
else else
xp=psi xp=psi
yp=phim yp=phim
xm=-psi xm=-psi
ym=phim ym=phim
x0=zero x0=0
y0=phim y0=phim
end if end if
call zetac (xp,yp,zrp,zip,iflag) call zetac (xp,yp,zrp,zip,iflag)
@ -1320,7 +1320,7 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error)
! !
cf12=czero cf12=czero
if (alpha.ge.0) then if (alpha.ge.0) then
if (alpha.ne.zero) cf12=-(czp+czm)/(2.0_wp_*phim) if (alpha.ne.0) cf12=-(czp+czm)/(2.0_wp_*phim)
else else
cf12=-im*(czp+czm)/(2.0_wp_*phim) cf12=-im*(czp+czm)/(2.0_wp_*phim)
end if end if

191
src/gray_core.f90
View File

@ -131,10 +131,10 @@ contains
allocate(results%jcd(params%output%nrho)) allocate(results%jcd(params%output%nrho))
! ...and initialise them ! ...and initialise them
results%pabs = zero results%pabs = 0
results%icd = zero results%icd = 0
results%dpdv = zero results%dpdv = 0
results%jcd = zero results%jcd = 0
! ========= set environment END ========= ! ========= set environment END =========
! Pre-determinted tables ! Pre-determinted tables
@ -180,15 +180,15 @@ contains
write (msg, '("pass: ",g0)') ip write (msg, '("pass: ",g0)') ip
call log_info(msg, mod='gray_core', proc='gray_main') call log_info(msg, mod='gray_core', proc='gray_main')
pabs_pass = zero pabs_pass = 0
icd_pass = zero icd_pass = 0
istop_pass = 0 ! stop flag for current pass istop_pass = 0 ! stop flag for current pass
nbeam_pass = 2*nbeam_pass ! max n of beams in current pass nbeam_pass = 2*nbeam_pass ! max n of beams in current pass
if(ip > 1) then if(ip > 1) then
du1 = zero du1 = 0
gri = zero gri = 0
ggri = zero ggri = 0
if(ip == params%raytracing%ipass) cpl = [zero, zero] ! no successive passes if(ip == params%raytracing%ipass) cpl = [zero, zero] ! no successive passes
end if end if
@ -215,15 +215,22 @@ contains
cycle cycle
end if end if
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) psjki = 0
ppabs = 0
ccci = 0
tau0 = 0
alphaabs0 = 0
dids0 = 0
ccci0 = 0
iiv = 1
if(ip == 1) then ! 1st pass if(ip == 1) then ! 1st pass
igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass
tau1 = zero ! * tau from previous passes tau1 = 0 ! * tau from previous passes
etau1 = one etau1 = 1
cpl1 = one ! * coupling from previous passes cpl1 = 1 ! * coupling from previous passes
lgcpl1 = zero lgcpl1 = 0
p0ray = p0jk ! * initial beam power p0ray = p0jk ! * initial beam power
call compute_initial_conds(params%raytracing, params%antenna, & ! * initial conditions call compute_initial_conds(params%raytracing, params%antenna, & ! * initial conditions
@ -469,10 +476,10 @@ contains
end if end if
end block end block
else else
tekev=zero tekev=0
alpha=zero alpha=0
didp=zero didp=0
anprim=zero anprim=0
anprre=anpr anprre=anpr
nharm=0 nharm=0
nhf=0 nhf=0
@ -568,15 +575,15 @@ contains
cpl_beam1 = sum(& cpl_beam1 = sum(&
p0ray * exp(-tau0) * cpls(:,child_index_rt)/cpl1, MASK=iop > 2) / & p0ray * exp(-tau0) * cpls(:,child_index_rt)/cpl1, MASK=iop > 2) / &
sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays) sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays)
cpl_beam2 = one - cpl_beam1 ! * average X-mode coupling for next beam cpl_beam2 = 1 - cpl_beam1 ! * average X-mode coupling for next beam
if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam
cpl_cbeam1 = cpls(1,child_index_rt)/cpl1(1) cpl_cbeam1 = cpls(1,child_index_rt)/cpl1(1)
cpl_cbeam2 = one - cpl_cbeam1 cpl_cbeam2 = 1 - cpl_cbeam1
end if end if
else ! last pass OR no ray re-entered plasma else ! last pass OR no ray re-entered plasma
cpl_beam1 = zero cpl_beam1 = 0
cpl_beam2 = zero cpl_beam2 = 0
end if end if
! print final results for pass on screen ! print final results for pass on screen
@ -661,36 +668,6 @@ contains
end subroutine gray_main end subroutine gray_main
subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
use const_and_precisions, only : zero
! arguments
real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci
real(wp_), dimension(:), intent(out) :: tau0,alphaabs0,dids0,ccci0
integer, dimension(:), intent(out) :: iiv
!! common/external functions/variables
! integer :: jclosest
! real(wp_), dimension(3) :: anwcl,xwcl
!
! common/refln/anwcl,xwcl,jclosest
!
! jclosest=nrayr+1
! anwcl(1:3)=0.0_wp_
! xwcl(1:3)=0.0_wp_
psjki = zero
ppabs = zero
ccci = zero
tau0 = zero
alphaabs0 = zero
dids0 = zero
ccci0 = zero
iiv = 1
end subroutine vectinit
subroutine compute_initial_conds(rtx, beam, N_c, k0, y, yp, dist, & subroutine compute_initial_conds(rtx, beam, N_c, k0, y, yp, dist, &
pos, grad_u, grad, hess) pos, grad_u, grad, hess)
! Computes the initial conditions for tracing a beam ! Computes the initial conditions for tracing a beam
@ -1007,7 +984,6 @@ contains
subroutine gradi_upd(params, ywrk, ak0, xc, du1, gri, ggri, error) subroutine gradi_upd(params, ywrk, ak0, xc, du1, gri, ggri, error)
use const_and_precisions, only : zero, half
use gray_params, only : raytracing_parameters use gray_params, only : raytracing_parameters
use gray_errors, only : gray_error, unstable_beam, raise_error use gray_errors, only : gray_error, unstable_beam, raise_error
@ -1060,7 +1036,7 @@ contains
call solg0(dxv1,dxv2,dxv3,dgu) call solg0(dxv1,dxv2,dxv3,dgu)
du1(:,k,j) = dgu du1(:,k,j) = dgu
end do end do
gri(:,1) = zero gri(:,1) = 0
! compute grad u1 and grad(S_I) for all the other rays ! compute grad u1 and grad(S_I) for all the other rays
if (params%nrayr > 1) then if (params%nrayr > 1) then
@ -1125,7 +1101,7 @@ contains
end do end do
! compute derivatives of grad u and grad(S_I) for rays jk>1 ! compute derivatives of grad u and grad(S_I) for rays jk>1
ggri(:,:,1) = zero ggri(:,:,1) = 0
jm=1 jm=1
j=2 j=2
k=0 k=0
@ -1160,9 +1136,9 @@ contains
uxx = dgg(1,1) uxx = dgg(1,1)
uyy = dgg(2,2) uyy = dgg(2,2)
uzz = dgg(3,3) uzz = dgg(3,3)
uxy = (dgg(1,2) + dgg(2,1))*half uxy = (dgg(1,2) + dgg(2,1))/2
uxz = (dgg(1,3) + dgg(3,1))*half uxz = (dgg(1,3) + dgg(3,1))/2
uyz = (dgg(2,3) + dgg(3,2))*half uyz = (dgg(2,3) + dgg(3,2))/2
! derivatives of S_I and Grad(S_I) ! derivatives of S_I and Grad(S_I)
gx = ux*dffiu gx = ux*dffiu
@ -1243,7 +1219,7 @@ contains
subroutine plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & subroutine plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, &
bv, derbv, xg, yg, derxg, deryg, psinv_opt) bv, derbv, xg, yg, derxg, deryg, psinv_opt)
use const_and_precisions, only : zero, cm use const_and_precisions, only : cm
use gray_equil, only : abstract_equil, vacuum use gray_equil, only : abstract_equil, vacuum
use gray_plasma, only : abstract_plasma use gray_plasma, only : abstract_plasma
@ -1267,15 +1243,15 @@ contains
real(wp_) :: brr,bphi,bzz,dxgdpsi real(wp_) :: brr,bphi,bzz,dxgdpsi
real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddensdpsi real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddensdpsi
xg = zero xg = 0
yg = 99._wp_ yg = 99._wp_
psinv = -1._wp_ psinv = -1._wp_
dens = zero dens = 0
btot = zero btot = 0
derxg = zero derxg = 0
deryg = zero deryg = 0
bv = zero bv = 0
derbv = zero derbv = 0
select type (equil) select type (equil)
type is (vacuum) type is (vacuum)
@ -1284,11 +1260,11 @@ contains
return return
end select end select
dbtot = zero dbtot = 0
dbv = zero dbv = 0
dbvcdc = zero dbvcdc = 0
dbvcdc = zero dbvcdc = 0
dbvdc = zero dbvdc = 0
xx = xv(1) xx = xv(1)
yy = xv(2) yy = xv(2)
@ -1315,7 +1291,7 @@ contains
if (present(psinv_opt)) psinv_opt = psinv if (present(psinv_opt)) psinv_opt = psinv
! compute yg and derivative ! compute yg and derivative
if(psinv < zero) then if(psinv < 0) then
bphi = fpolv/rrm bphi = fpolv/rrm
btot = abs(bphi) btot = abs(bphi)
yg = btot/bres yg = btot/bres
@ -1405,7 +1381,7 @@ contains
! (`ywppla_upd` suborutine); while the optional ones are used for ! (`ywppla_upd` suborutine); while the optional ones are used for
! computing the absoprtion and current drive. ! computing the absoprtion and current drive.
use const_and_precisions, only : zero, one, half, two use const_and_precisions, only : zero, half
use gray_params, only : gray_parameters, STEP_ARCLEN, & use gray_params, only : gray_parameters, STEP_ARCLEN, &
STEP_TIME, STEP_PHASE STEP_TIME, STEP_PHASE
@ -1461,24 +1437,24 @@ contains
! Shorthands used in the expressions below ! Shorthands used in the expressions below
yg2 = yg**2 ! Y² yg2 = yg**2 ! Y²
anpl2 = anpl**2 ! N∥² anpl2 = anpl**2 ! N∥²
dnl = one - anpl2 ! 1 - N∥² dnl = 1 - anpl2 ! 1 - N∥²
duh = one - xg - yg2 ! UH denom (duh=0 on the upper-hybrid resonance) duh = 1 - xg - yg2 ! UH denom (duh=0 on the upper-hybrid resonance)
! Compute/copy optional outputs ! Compute/copy optional outputs
if (present(anpr)) anpr = sqrt(max(an2 - anpl2, zero)) ! N⊥ if (present(anpr)) anpr = sqrt(max(an2 - anpl2, zero)) ! N⊥
if (present(anpl_)) anpl_ = anpl ! N∥ if (present(anpl_)) anpl_ = anpl ! N∥
an2s = one an2s = 1
dan2sdxg = zero dan2sdxg = 0
dan2sdyg = zero dan2sdyg = 0
dan2sdnpl = zero dan2sdnpl = 0
del = zero del = 0
fdia = zero fdia = 0
dfdiadnpl = zero dfdiadnpl = 0
dfdiadxg = zero dfdiadxg = 0
dfdiadyg = zero dfdiadyg = 0
if(xg > zero) then if(xg > 0) then
! Derivatives of the cold plasma refractive index ! Derivatives of the cold plasma refractive index
! !
! N²s = 1 - X - XY²⋅(1 + N∥² ± √Δ)/[2(1 - X - Y²)] ! N²s = 1 - X - XY²⋅(1 + N∥² ± √Δ)/[2(1 - X - Y²)]
@ -1487,19 +1463,19 @@ contains
! + for the X mode, - for the O mode ! + for the X mode, - for the O mode
! √Δ ! √Δ
del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2) del = sqrt(dnl**2 + 4.0_wp_*anpl2*(1 - xg)/yg2)
! ∂(N²s)/∂X ! ∂(N²s)/∂X
! Note: this term is nonzero for X=0, but it multiplies terms ! Note: this term is nonzero for X=0, but it multiplies terms
! proportional to X or ∂X/∂ψ which are zero outside the plasma. ! proportional to X or ∂X/∂ψ which are zero outside the plasma.
dan2sdxg = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2 & dan2sdxg = - half*yg2*(1 - yg2)*(1 + anpl2 + sox*del)/duh**2 &
+ sox*xg*anpl2/(del*duh) - one + sox*xg*anpl2/(del*duh) - 1
! ∂(N²s)/∂Y ! ∂(N²s)/∂Y
dan2sdyg = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2 & dan2sdyg = - xg*yg*(1 - xg)*(1 + anpl2 + sox*del)/duh**2 &
+ two*sox*xg*(one - xg)*anpl2/(yg*del*duh) + 2*sox*xg*(1 - xg)*anpl2/(yg*del*duh)
! ∂(N²s)/∂N∥ ! ∂(N²s)/∂N∥
dan2sdnpl = - xg*yg2*anpl/duh & dan2sdnpl = - xg*yg2*anpl/duh &
- sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh) - sox*xg*anpl*(2*(1 - xg) - yg2*dnl)/(del*duh)
if(igrad > 0) then if(igrad > 0) then
! Derivatives used in the complex eikonal terms (beamtracing only) ! Derivatives used in the complex eikonal terms (beamtracing only)
@ -1507,28 +1483,27 @@ contains
real(wp_) :: ddelnpl2, ddelnpl2x, ddelnpl2y, derdel real(wp_) :: ddelnpl2, ddelnpl2x, ddelnpl2y, derdel
! ∂²Δ/∂N∥² ! ∂²Δ/∂N∥²
ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) & ddelnpl2 = 2*(2*(1 - xg)*(1 + 3.0_wp_*anpl2**2) &
- yg2*dnl**3)/yg2/del**3 - yg2*dnl**3)/yg2/del**3
! ∂²(N²s)/∂N∥² ! ∂²(N²s)/∂N∥²
fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh fdia = - xg*yg2*(1 + half*sox*ddelnpl2)/duh
! Intermediates results used right below ! Intermediates results used right below
derdel = two*(one - xg)*anpl2*(one + 3.0_wp_*anpl2**2) & derdel = 2*(1 - xg)*anpl2*(1 + 3*anpl2**2) &
- dnl**2*(one + 3.0_wp_*anpl2)*yg2 - dnl**2*(1 + 3*anpl2)*yg2
derdel = 4.0_wp_*derdel/(yg*del)**5 derdel = 4*derdel/(yg*del)**5
ddelnpl2y = two*(one - xg)*derdel ddelnpl2y = 2*(1 - xg)*derdel
ddelnpl2x = yg*derdel ddelnpl2x = yg*derdel
! ∂³(N²s)/∂N∥³ ! ∂³(N²s)/∂N∥³
dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2) & dfdiadnpl = 24*sox*xg*(1 - xg)*anpl*(1 - anpl2**2)/(yg2*del**5)
/(yg2*del**5)
! ∂³(N²s)/∂N∥²∂X ! ∂³(N²s)/∂N∥²∂X
dfdiadxg = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2) & dfdiadxg = - yg2*(1 - yg2)/duh**2 - sox*yg2*((1 - yg2) &
*ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2) *ddelnpl2 + xg*duh*ddelnpl2x)/(2*duh**2)
! ∂³(N²s)/∂N∥²∂Y ! ∂³(N²s)/∂N∥²∂Y
dfdiadyg = - two*yg*xg*(one - xg)/duh**2 & dfdiadyg = - 2*yg*xg*(1 - xg)/duh**2 &
- sox*xg*yg*(two*(one - xg)*ddelnpl2 & - sox*xg*yg*(2*(1 - xg)*ddelnpl2 &
+ yg*duh*ddelnpl2y)/(two*duh**2) + yg*duh*ddelnpl2y)/(2*duh**2)
end block end block
end if end if
end if end if
@ -1541,7 +1516,7 @@ contains
! ∂Λ/∂N̅ = 2N̅ - ∂(N²s)/∂N̅ = 2N̅ - ∂(N²s)/∂N∥ b̅ ! ∂Λ/∂N̅ = 2N̅ - ∂(N²s)/∂N̅ = 2N̅ - ∂(N²s)/∂N∥ b̅
! Note: we used the identity ∇f(v̅⋅b̅) = f' ∇(v̅⋅b̅) = f'b̅. ! Note: we used the identity ∇f(v̅⋅b̅) = f' ∇(v̅⋅b̅) = f'b̅.
derdnv = two*anv - dan2sdnpl*bv derdnv = 2*anv - dan2sdnpl*bv
! ∂Λ/∂ω = ∂N²/∂ω - ∂N²s/∂X⋅∂X/∂ω - ∂N²s/∂Y⋅∂Y/∂ω - ∂N²s/∂N∥⋅∂N∥/∂ω ! ∂Λ/∂ω = ∂N²/∂ω - ∂N²s/∂X⋅∂X/∂ω - ∂N²s/∂Y⋅∂Y/∂ω - ∂N²s/∂N∥⋅∂N∥/∂ω
! Notes: 1. N depends on ω: N²=c²k²/ω² ⇒ ∂N²/∂ω = -2N²/ω ! Notes: 1. N depends on ω: N²=c²k²/ω² ⇒ ∂N²/∂ω = -2N²/ω
@ -1549,7 +1524,7 @@ contains
! 2. derdom is actually ω⋅∂Λ/∂ω, see below for the reason. ! 2. derdom is actually ω⋅∂Λ/∂ω, see below for the reason.
! 3. N gains a dependency on ω because Λ(∇S, ω) is computed ! 3. N gains a dependency on ω because Λ(∇S, ω) is computed
! on the constrains Λ=0. ! on the constrains Λ=0.
derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl derdom = -2*an2 + 2*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl
if (igrad > 0) then if (igrad > 0) then
! Complex eikonal terms added to the above expressions ! Complex eikonal terms added to the above expressions
@ -1576,7 +1551,7 @@ contains
! ∂Λ/∂ω += ∂|∇S_I|²/∂ω + ½∂(b⋅∇S_I)²/∂ω + ½(b̅⋅∇S_I)² ∂/∂ω (∂²N²s/∂N∥²) ! ∂Λ/∂ω += ∂|∇S_I|²/∂ω + ½∂(b⋅∇S_I)²/∂ω + ½(b̅⋅∇S_I)² ∂/∂ω (∂²N²s/∂N∥²)
! Note: as above ∇S_I gains a dependency on ω ! Note: as above ∇S_I gains a dependency on ω
derdom = derdom + two*gr2 - bdotgr**2 & derdom = derdom + 2*gr2 - bdotgr**2 &
* (fdia + xg*dfdiadxg + half*yg*dfdiadyg & * (fdia + xg*dfdiadxg + half*yg*dfdiadyg &
+ half*anpl*dfdiadnpl) + half*anpl*dfdiadnpl)
end block end block
@ -1651,7 +1626,7 @@ contains
if (present(ddr) .or. present(ddi)) then if (present(ddr) .or. present(ddi)) then
! Dispersion relation (geometric optics) ! Dispersion relation (geometric optics)
! ddr ← Λ = N² - N²s(X,Y,N∥) = 0 ! ddr ← Λ = N² - N²s(X,Y,N∥) = 0
an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh an2s = 1 - xg - half*xg*yg2*(1 + anpl2 + sox*del)/duh
ddr = an2 - an2s ddr = an2 - an2s
end if end if

2
src/gray_plasma.f90
View File

@ -185,7 +185,7 @@ contains
ddens = self%dens_spline%deriv(psin) ddens = self%dens_spline%deriv(psin)
! Evaluate the spline 1st derivative ! Evaluate the spline 1st derivative
if (abs(dens) < 1.0e-10_wp_) dens = zero if (abs(dens) < 1.0e-10_wp_) dens = 0
else else
! Use a C² polynomial extension outside (ψ > ψ₀) ! Use a C² polynomial extension outside (ψ > ψ₀)

22
src/reflections.f90
View File

@ -1,5 +1,5 @@
module reflections module reflections
use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge, zero, one use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge
implicit none implicit none
@ -16,7 +16,7 @@ subroutine reflect(ki,nsurf,ko)
real(wp_), intent(out), dimension(3) :: ko real(wp_), intent(out), dimension(3) :: ko
real(wp_) :: twokn,norm2 real(wp_) :: twokn,norm2
norm2 = dot_product(nsurf,nsurf) norm2 = dot_product(nsurf,nsurf)
if (norm2>zero) then if (norm2>0) then
twokn = 2.0_wp_*dot_product(ki,nsurf)/norm2 twokn = 2.0_wp_*dot_product(ki,nsurf)/norm2
ko=ki-twokn*nsurf ko=ki-twokn*nsurf
else else
@ -42,20 +42,20 @@ subroutine inters_linewall(xv, kv, wall, sint, normw)
sint=comp_huge sint=comp_huge
iint=0 iint=0
normw=zero normw=0
do i=1, size(wall%R)-1 do i=1, size(wall%R)-1
! search intersections with i-th wall segment ! search intersections with i-th wall segment
call linecone_coord(xv, kv, wall%R(i:i+1), wall%z(i:i+1), si, ti, ni) call linecone_coord(xv, kv, wall%R(i:i+1), wall%z(i:i+1), si, ti, ni)
! discard solutions with s<=0 ! discard solutions with s<=0
first=ni+1 first=ni+1
do j=1,ni do j=1,ni
if (si(j)>zero) then if (si(j)>0) then
first=j first=j
exit exit
end if end if
end do end do
do j=first,ni do j=first,ni
if ((si(j)<sint .or. iint==0) .and. ti(j)>=zero .and. ti(j)<=one) then if ((si(j)<sint .or. iint==0) .and. ti(j)>=0 .and. ti(j)<=1) then
! check intersection is in r,z range and keep the closest ! check intersection is in r,z range and keep the closest
sint = si(j) sint = si(j)
iint = i iint = i
@ -72,7 +72,7 @@ subroutine inters_linewall(xv, kv, wall, sint, normw)
l = hypot(drw, dzw) l = hypot(drw, dzw)
kxy = norm2(kv(1:2)) kxy = norm2(kv(1:2))
normw(3) = -drw/l normw(3) = -drw/l
if (rint>zero) then if (rint>0) then
normw(1) = xint/rint*dzw/l normw(1) = xint/rint*dzw/l
normw(2) = yint/rint*dzw/l normw(2) = yint/rint*dzw/l
else else
@ -80,7 +80,7 @@ subroutine inters_linewall(xv, kv, wall, sint, normw)
normw(2) = kv(2)/kxy*dzw/l normw(2) = kv(2)/kxy*dzw/l
end if end if
! reverse normal if k.n>0 ! reverse normal if k.n>0
if (dot_product(normw,kv)>zero) normw=-normw if (dot_product(normw,kv)>0) normw=-normw
end subroutine inters_linewall end subroutine inters_linewall
@ -168,8 +168,8 @@ subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr)
dyb = yb(2)-yb(1) dyb = yb(2)-yb(1)
crossprod = dxb*dya - dxa*dyb crossprod = dxb*dya - dxa*dyb
if (abs(crossprod)<comp_tiny) then if (abs(crossprod)<comp_tiny) then
s = zero s = 0
t = zero t = 0
ierr = 1 ierr = 1
else else
s = (dyb*(xa(1)-xb(1)) - dxb*(ya(1)-yb(1)))/crossprod s = (dyb*(xa(1)-xb(1)) - dxb*(ya(1)-yb(1)))/crossprod
@ -187,8 +187,8 @@ function interssegm(xa,ya,xb,yb)
integer :: ierr integer :: ierr
interssegm = .false. interssegm = .false.
call interssegm_coord(xa,ya,xb,yb,s,t,ierr) call interssegm_coord(xa,ya,xb,yb,s,t,ierr)
if (ierr==0 .and. s>=zero .and. s<=one .and. & if (ierr==0 .and. s>=0 .and. s<=1 .and. &
t>=zero .and. t<=one) interssegm = .true. t>=0 .and. t<=1) interssegm = .true.
end function interssegm end function interssegm

26
src/vendor/eierf.f90 vendored
View File

@ -180,7 +180,7 @@ contains
1.99999999999048104167_wp_/) 1.99999999999048104167_wp_/)
!---------------------------------------------------------------------- !----------------------------------------------------------------------
x = arg x = arg
if (x == zero) then if (x == 0) then
ei = -xinf ei = -xinf
if (intt == 2) ei = -ei if (intt == 2) ei = -ei
else if ((x < zero) .or. (intt == 2)) then else if ((x < zero) .or. (intt == 2)) then
@ -209,7 +209,7 @@ contains
if (intt /= 3) ei = ei * exp(-y) if (intt /= 3) ei = ei * exp(-y)
else else
if ((y > xbig) .and. (intt < 3)) then if ((y > xbig) .and. (intt < 3)) then
ei = zero ei = 0
else else
w = one / y w = one / y
sump = e(1) sump = e(1)
@ -231,8 +231,8 @@ contains
!---------------------------------------------------------------------- !----------------------------------------------------------------------
t = x + x t = x + x
t = t / three - two t = t / three - two
px(1) = zero px(1) = 0
qx(1) = zero qx(1) = 0
px(2) = p(1) px(2) = p(1)
qx(2) = q(1) qx(2) = q(1)
do i = 2, 9 do i = 2, 9
@ -262,14 +262,14 @@ contains
if (intt == 3) ei = exp(-x) * ei if (intt == 3) ei = exp(-x) * ei
end if end if
else if (x < twelve) then else if (x < twelve) then
frac = zero frac = 0
do i = 1, 9 do i = 1, 9
frac = s(i) / (r(i) + x + frac) frac = s(i) / (r(i) + x + frac)
end do end do
ei = (r(10) + frac) / x ei = (r(10) + frac) / x
if (intt /= 3) ei = ei * exp(x) if (intt /= 3) ei = ei * exp(x)
else if (x <= two4) then else if (x <= two4) then
frac = zero frac = 0
do i = 1, 9 do i = 1, 9
frac = q1(i) / (p1(i) + x + frac) frac = q1(i) / (p1(i) + x + frac)
end do end do
@ -280,7 +280,7 @@ contains
ei = xinf ei = xinf
else else
y = one / x y = one / x
frac = zero frac = 0
do i = 1, 9 do i = 1, 9
frac = q2(i) / (p2(i) + x + frac) frac = q2(i) / (p2(i) + x + frac)
end do end do
@ -532,7 +532,7 @@ contains
1.99999999999048104167_wp_/) 1.99999999999048104167_wp_/)
!---------------------------------------------------------------------- !----------------------------------------------------------------------
x = arg x = arg
if (x == zero) then if (x == 0) then
ei = -xinf ei = -xinf
else if ((x < zero)) then else if ((x < zero)) then
!---------------------------------------------------------------------- !----------------------------------------------------------------------
@ -574,8 +574,8 @@ contains
!---------------------------------------------------------------------- !----------------------------------------------------------------------
t = x + x t = x + x
t = t / three - two t = t / three - two
px(1) = zero px(1) = 0
qx(1) = zero qx(1) = 0
px(2) = p(1) px(2) = p(1)
qx(2) = q(1) qx(2) = q(1)
do i = 2, 9 do i = 2, 9
@ -603,20 +603,20 @@ contains
ei = exp(-x) * (sump / (sumq*(x+x0)) + frac) * xmx0 ei = exp(-x) * (sump / (sumq*(x+x0)) + frac) * xmx0
end if end if
else if (x < twelve) then else if (x < twelve) then
frac = zero frac = 0
do i = 1, 9 do i = 1, 9
frac = s(i) / (r(i) + x + frac) frac = s(i) / (r(i) + x + frac)
end do end do
ei = (r(10) + frac) / x ei = (r(10) + frac) / x
else if (x <= two4) then else if (x <= two4) then
frac = zero frac = 0
do i = 1, 9 do i = 1, 9
frac = q1(i) / (p1(i) + x + frac) frac = q1(i) / (p1(i) + x + frac)
end do end do
ei = (p1(10) + frac) / x ei = (p1(10) + frac) / x
else else
y = one / x y = one / x
frac = zero frac = 0
do i = 1, 9 do i = 1, 9
frac = q2(i) / (p2(i) + x + frac) frac = q2(i) / (p2(i) + x + frac)
end do end do

156
src/vendor/minpack.f90 vendored
View File

@ -48,7 +48,6 @@ module minpack
contains contains
pure subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa) pure subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa)
use const_and_precisions, only : zero, one
! arguments ! arguments
integer, intent(in) :: n, ldfjac, lwa integer, intent(in) :: n, ldfjac, lwa
integer, intent(out) :: info integer, intent(out) :: info
@ -166,7 +165,7 @@ contains
! !
! check the input parameters for errors. ! check the input parameters for errors.
! !
if (n <= 0 .or. ldfjac < n .or. tol < zero & if (n <= 0 .or. ldfjac < n .or. tol < 0 &
.or. lwa < (n*(n + 13))/2) return .or. lwa < (n*(n + 13))/2) return
! !
! call hybrj. ! call hybrj.
@ -175,7 +174,7 @@ contains
xtol = tol xtol = tol
mode = 2 mode = 2
do j = 1, n do j = 1, n
wa(j) = one wa(j) = 1
end do end do
nprint = 0 nprint = 0
lr = (n*(n + 1))/2 lr = (n*(n + 1))/2
@ -188,7 +187,7 @@ contains
pure subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, & pure subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, &
factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, & factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, &
wa3,wa4) wa3,wa4)
use const_and_precisions, only : zero, one, epsmch=>comp_eps use const_and_precisions, only : epsmch=>comp_eps
use, intrinsic :: ieee_exceptions, only : ieee_get_halting_mode, & use, intrinsic :: ieee_exceptions, only : ieee_get_halting_mode, &
ieee_set_halting_mode, & ieee_set_halting_mode, &
ieee_invalid ieee_invalid
@ -373,12 +372,12 @@ contains
! !
! check the input parameters for errors. ! check the input parameters for errors.
! !
if (n <= 0 .or. ldfjac < n .or. xtol < zero & if (n <= 0 .or. ldfjac < n .or. xtol < 0 &
.or. maxfev <= 0 .or. factor <= zero & .or. maxfev <= 0 .or. factor <= 0 &
.or. lr < (n*(n + 1))/2) go to 300 .or. lr < (n*(n + 1))/2) go to 300
if (mode == 2) then if (mode == 2) then
do j = 1, n do j = 1, n
if (diag(j) <= zero) go to 300 if (diag(j) <= 0) go to 300
end do end do
end if end if
! !
@ -422,7 +421,7 @@ contains
if (mode /= 2) then if (mode /= 2) then
do j = 1, n do j = 1, n
diag(j) = wa2(j) diag(j) = wa2(j)
if (wa2(j) == zero) diag(j) = one if (wa2(j) == 0) diag(j) = 1
end do end do
end if end if
! !
@ -434,7 +433,7 @@ contains
end do end do
xnorm = enorm(n,wa3) xnorm = enorm(n,wa3)
delta = factor*xnorm delta = factor*xnorm
if (delta == zero) delta = factor if (delta == 0) delta = factor
end if end if
! !
! form (q transpose)*fvec and store in qtf. ! form (q transpose)*fvec and store in qtf.
@ -443,8 +442,8 @@ contains
qtf(i) = fvec(i) qtf(i) = fvec(i)
end do end do
do j = 1, n do j = 1, n
if (fjac(j,j) /= zero) then if (fjac(j,j) /= 0) then
summ = zero summ = 0
do i = j, n do i = j, n
summ = summ + fjac(i,j)*qtf(i) summ = summ + fjac(i,j)*qtf(i)
end do end do
@ -466,7 +465,7 @@ contains
l = l + n - i l = l + n - i
end do end do
r(l) = wa1(j) r(l) = wa1(j)
if (wa1(j) == zero) sing = .true. if (wa1(j) == 0) sing = .true.
end do end do
! !
! accumulate the orthogonal factor in fjac. ! accumulate the orthogonal factor in fjac.
@ -520,14 +519,14 @@ contains
! !
! compute the scaled actual reduction. ! compute the scaled actual reduction.
! !
actred = -one actred = -1
if (fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2 if (fnorm1 < fnorm) actred = 1 - (fnorm1/fnorm)**2
! !
! compute the scaled predicted reduction. ! compute the scaled predicted reduction.
! !
l = 1 l = 1
do i = 1, n do i = 1, n
summ = zero summ = 0
do j = i, n do j = i, n
summ = summ + r(l)*wa1(j) summ = summ + r(l)*wa1(j)
l = l + 1 l = l + 1
@ -535,14 +534,14 @@ contains
wa3(i) = qtf(i) + summ wa3(i) = qtf(i) + summ
end do end do
temp = enorm(n,wa3) temp = enorm(n,wa3)
prered = zero prered = 0
if (temp < fnorm) prered = one - (temp/fnorm)**2 if (temp < fnorm) prered = 1 - (temp/fnorm)**2
! !
! compute the ratio of the actual to the predicted ! compute the ratio of the actual to the predicted
! reduction. ! reduction.
! !
ratio = zero ratio = 0
if (prered > zero) ratio = actred/prered if (prered > 0) ratio = actred/prered
! !
! update the step bound. ! update the step bound.
! !
@ -554,7 +553,7 @@ contains
ncfail = 0 ncfail = 0
ncsuc = ncsuc + 1 ncsuc = ncsuc + 1
if (ratio >= p5 .or. ncsuc > 1) delta = dmax1(delta,pnorm/p5) if (ratio >= p5 .or. ncsuc > 1) delta = dmax1(delta,pnorm/p5)
if (abs(ratio-one) <= p1) delta = pnorm/p5 if (abs(ratio-1) <= p1) delta = pnorm/p5
end if end if
! !
! test for successful iteration. ! test for successful iteration.
@ -582,7 +581,7 @@ contains
! !
! test for convergence. ! test for convergence.
! !
if (delta <= xtol*xnorm .or. fnorm == zero) info = 1 if (delta <= xtol*xnorm .or. fnorm == 0) info = 1
if (info /= 0) go to 300 if (info /= 0) go to 300
! !
! tests for termination and stringent tolerances. ! tests for termination and stringent tolerances.
@ -601,7 +600,7 @@ contains
! and update qtf if necessary. ! and update qtf if necessary.
! !
do j = 1, n do j = 1, n
summ = zero summ = 0
do i = 1, n do i = 1, n
summ = summ + fjac(i,j)*wa4(i) summ = summ + fjac(i,j)*wa4(i)
end do end do
@ -638,7 +637,7 @@ contains
end subroutine hybrj end subroutine hybrj
pure subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2) pure subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
use const_and_precisions, only : zero, one, epsmch=>comp_eps use const_and_precisions, only : epsmch=>comp_eps
! arguments ! arguments
integer, intent(in) :: n, lr integer, intent(in) :: n, lr
real(wp_), intent(in) :: delta, r(lr), diag(n), qtb(n) real(wp_), intent(in) :: delta, r(lr), diag(n), qtb(n)
@ -712,20 +711,20 @@ contains
jp1 = j + 1 jp1 = j + 1
jj = jj - k jj = jj - k
l = jj + 1 l = jj + 1
summ = zero summ = 0
do i = jp1, n do i = jp1, n
summ = summ + r(l)*x(i) summ = summ + r(l)*x(i)
l = l + 1 l = l + 1
end do end do
temp = r(jj) temp = r(jj)
if (temp == zero) then if (temp == 0) then
l = j l = j
do i = 1, j do i = 1, j
temp = dmax1(temp,abs(r(l))) temp = dmax1(temp,abs(r(l)))
l = l + n - i l = l + n - i
end do end do
temp = epsmch*temp temp = epsmch*temp
if (temp == zero) temp = epsmch if (temp == 0) temp = epsmch
end if end if
x(j) = (qtb(j) - summ)/temp x(j) = (qtb(j) - summ)/temp
end do end do
@ -733,7 +732,7 @@ contains
! test whether the gauss-newton direction is acceptable. ! test whether the gauss-newton direction is acceptable.
! !
do j = 1, n do j = 1, n
wa1(j) = zero wa1(j) = 0
wa2(j) = diag(j)*x(j) wa2(j) = diag(j)*x(j)
end do end do
qnorm = enorm(n,wa2) qnorm = enorm(n,wa2)
@ -756,9 +755,9 @@ contains
! the special case in which the scaled gradient is zero. ! the special case in which the scaled gradient is zero.
! !
gnorm = enorm(n,wa1) gnorm = enorm(n,wa1)
sgnorm = zero sgnorm = 0
alpha = delta/qnorm alpha = delta/qnorm
if (gnorm /= zero) then if (gnorm /= 0) then
! !
! calculate the point along the scaled gradient ! calculate the point along the scaled gradient
! at which the quadratic is minimized. ! at which the quadratic is minimized.
@ -768,7 +767,7 @@ contains
end do end do
l = 1 l = 1
do j = 1, n do j = 1, n
summ = zero summ = 0
do i = j, n do i = j, n
summ = summ + r(l)*wa1(i) summ = summ + r(l)*wa1(i)
l = l + 1 l = l + 1
@ -780,7 +779,7 @@ contains
! !
! test whether the scaled gradient direction is acceptable. ! test whether the scaled gradient direction is acceptable.
! !
alpha = zero alpha = 0
if (sgnorm < delta) then if (sgnorm < delta) then
! !
! the scaled gradient direction is not acceptable. ! the scaled gradient direction is not acceptable.
@ -791,22 +790,21 @@ contains
temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta) temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
temp = temp - (delta/qnorm)*(sgnorm/delta)**2 & temp = temp - (delta/qnorm)*(sgnorm/delta)**2 &
+ sqrt((temp-(delta/qnorm))**2 & + sqrt((temp-(delta/qnorm))**2 &
+(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2)) +(1-(delta/qnorm)**2)*(1-(sgnorm/delta)**2))
alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp alpha = ((delta/qnorm)*(1 - (sgnorm/delta)**2))/temp
end if end if
end if end if
! !
! form appropriate convex combination of the gauss-newton ! form appropriate convex combination of the gauss-newton
! direction and the scaled gradient direction. ! direction and the scaled gradient direction.
! !
temp = (one - alpha)*dmin1(sgnorm,delta) temp = (1 - alpha)*dmin1(sgnorm,delta)
do j = 1, n do j = 1, n
x(j) = temp*wa1(j) + alpha*x(j) x(j) = temp*wa1(j) + alpha*x(j)
end do end do
end subroutine dogleg end subroutine dogleg
pure function enorm(n,x) pure function enorm(n,x)
use const_and_precisions, only : zero, one
real(wp_) :: enorm real(wp_) :: enorm
integer, intent(in) :: n integer, intent(in) :: n
real(wp_), dimension(n), intent(in) :: x real(wp_), dimension(n), intent(in) :: x
@ -850,11 +848,11 @@ contains
integer :: i integer :: i
real(wp_) :: agiant,floatn,s1,s2,s3,xabs,x1max,x3max real(wp_) :: agiant,floatn,s1,s2,s3,xabs,x1max,x3max
real(wp_), parameter :: rdwarf=3.834e-20_wp_,rgiant=1.304e19_wp_ real(wp_), parameter :: rdwarf=3.834e-20_wp_,rgiant=1.304e19_wp_
s1 = zero s1 = 0
s2 = zero s2 = 0
s3 = zero s3 = 0
x1max = zero x1max = 0
x3max = zero x3max = 0
floatn = n floatn = n
agiant = rgiant/floatn agiant = rgiant/floatn
do i = 1, n do i = 1, n
@ -865,7 +863,7 @@ contains
! sum for large components. ! sum for large components.
! !
if (xabs > x1max) then if (xabs > x1max) then
s1 = one + s1*(x1max/xabs)**2 s1 = 1 + s1*(x1max/xabs)**2
x1max = xabs x1max = xabs
else else
s1 = s1 + (xabs/x1max)**2 s1 = s1 + (xabs/x1max)**2
@ -875,10 +873,10 @@ contains
! sum for small components. ! sum for small components.
! !
if (xabs > x3max) then if (xabs > x3max) then
s3 = one + s3*(x3max/xabs)**2 s3 = 1 + s3*(x3max/xabs)**2
x3max = xabs x3max = xabs
else else
if (xabs /= zero) s3 = s3 + (xabs/x3max)**2 if (xabs /= 0) s3 = s3 + (xabs/x3max)**2
end if end if
end if end if
else else
@ -891,11 +889,11 @@ contains
! !
! calculation of norm. ! calculation of norm.
! !
if (s1 /= zero) then if (s1 /= 0) then
enorm = x1max*sqrt(s1+(s2/x1max)/x1max) enorm = x1max*sqrt(s1+(s2/x1max)/x1max)
else else
if (s2 /= zero) then if (s2 /= 0) then
if (s2 >= x3max) enorm = sqrt(s2*(one+(x3max/s2)*(x3max*s3))) if (s2 >= x3max) enorm = sqrt(s2*(1+(x3max/s2)*(x3max*s3)))
if (s2 < x3max) enorm = sqrt(x3max*((s2/x3max)+(x3max*s3))) if (s2 < x3max) enorm = sqrt(x3max*((s2/x3max)+(x3max*s3)))
else else
enorm = x3max*sqrt(s3) enorm = x3max*sqrt(s3)
@ -904,7 +902,6 @@ contains
end function enorm end function enorm
pure subroutine qform(m,n,q,ldq,wa) pure subroutine qform(m,n,q,ldq,wa)
use const_and_precisions, only : zero, one
! arguments ! arguments
integer, intent(in) :: m,n,ldq integer, intent(in) :: m,n,ldq
real(wp_), intent(out) :: wa(m) real(wp_), intent(out) :: wa(m)
@ -956,7 +953,7 @@ contains
do j = 2, minmn do j = 2, minmn
jm1 = j - 1 jm1 = j - 1
do i = 1, jm1 do i = 1, jm1
q(i,j) = zero q(i,j) = 0
end do end do
end do end do
! !
@ -965,9 +962,9 @@ contains
np1 = n + 1 np1 = n + 1
do j = np1, m do j = np1, m
do i = 1, m do i = 1, m
q(i,j) = zero q(i,j) = 0
end do end do
q(j,j) = one q(j,j) = 1
end do end do
! !
! accumulate q from its factored form. ! accumulate q from its factored form.
@ -976,12 +973,12 @@ contains
k = minmn - l + 1 k = minmn - l + 1
do i = k, m do i = k, m
wa(i) = q(i,k) wa(i) = q(i,k)
q(i,k) = zero q(i,k) = 0
end do end do
q(k,k) = one q(k,k) = 1
if (wa(k) /= zero) then if (wa(k) /= 0) then
do j = k, m do j = k, m
summ = zero summ = 0
do i = k, m do i = k, m
summ = summ + q(i,j)*wa(i) summ = summ + q(i,j)*wa(i)
end do end do
@ -995,7 +992,7 @@ contains
end subroutine qform end subroutine qform
pure subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa) pure subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
use const_and_precisions, only : zero, one, epsmch=>comp_eps use const_and_precisions, only : zero, epsmch=>comp_eps
! arguments ! arguments
integer, intent(in) :: m, n, lda, lipvt integer, intent(in) :: m, n, lda, lipvt
integer, intent(out) :: ipvt(lipvt) integer, intent(out) :: ipvt(lipvt)
@ -1122,19 +1119,19 @@ contains
! j-th column of a to a multiple of the j-th unit vector. ! j-th column of a to a multiple of the j-th unit vector.
! !
ajnorm = enorm(m-j+1,a(j,j)) ajnorm = enorm(m-j+1,a(j,j))
if (ajnorm /= zero) then if (ajnorm /= 0) then
if (a(j,j) < zero) ajnorm = -ajnorm if (a(j,j) < 0) ajnorm = -ajnorm
do i = j, m do i = j, m
a(i,j) = a(i,j)/ajnorm a(i,j) = a(i,j)/ajnorm
end do end do
a(j,j) = a(j,j) + one a(j,j) = a(j,j) + 1
! !
! apply the transformation to the remaining columns ! apply the transformation to the remaining columns
! and update the norms. ! and update the norms.
! !
jp1 = j + 1 jp1 = j + 1
do k = jp1, n do k = jp1, n
summ = zero summ = 0
do i = j, m do i = j, m
summ = summ + a(i,j)*a(i,k) summ = summ + a(i,j)*a(i,k)
end do end do
@ -1142,9 +1139,9 @@ contains
do i = j, m do i = j, m
a(i,k) = a(i,k) - temp*a(i,j) a(i,k) = a(i,k) - temp*a(i,j)
end do end do
if (pivot .and. rdiag(k) /= zero) then if (pivot .and. rdiag(k) /= 0) then
temp = a(j,k)/rdiag(k) temp = a(j,k)/rdiag(k)
rdiag(k) = rdiag(k)*sqrt(dmax1(zero,one-temp**2)) rdiag(k) = rdiag(k)*sqrt(dmax1(zero, 1-temp**2))
if (p05*(rdiag(k)/wa(k))**2 <= epsmch) then if (p05*(rdiag(k)/wa(k))**2 <= epsmch) then
rdiag(k) = enorm(m-j,a(jp1,k)) rdiag(k) = enorm(m-j,a(jp1,k))
wa(k) = rdiag(k) wa(k) = rdiag(k)
@ -1157,7 +1154,6 @@ contains
end subroutine qrfac end subroutine qrfac
pure subroutine r1mpyq(m,n,a,lda,v,w) pure subroutine r1mpyq(m,n,a,lda,v,w)
use const_and_precisions, only : one
! arguments ! arguments
integer, intent(in) :: m, n, lda integer, intent(in) :: m, n, lda
real(wp_), intent(in) :: v(n),w(n) real(wp_), intent(in) :: v(n),w(n)
@ -1221,10 +1217,10 @@ contains
if (nm1 < 1) return if (nm1 < 1) return
do nmj = 1, nm1 do nmj = 1, nm1
j = n - nmj j = n - nmj
if (abs(v(j)) > one) cs = one/v(j) if (abs(v(j)) > 1) cs = 1/v(j)
if (abs(v(j)) > one) sn = sqrt(one-cs**2) if (abs(v(j)) > 1) sn = sqrt(1-cs**2)
if (abs(v(j)) <= one) sn = v(j) if (abs(v(j)) <= 1) sn = v(j)
if (abs(v(j)) <= one) cs = sqrt(one-sn**2) if (abs(v(j)) <= 1) cs = sqrt(1-sn**2)
do i = 1, m do i = 1, m
temp = cs*a(i,j) - sn*a(i,n) temp = cs*a(i,j) - sn*a(i,n)
a(i,n) = sn*a(i,j) + cs*a(i,n) a(i,n) = sn*a(i,j) + cs*a(i,n)
@ -1235,10 +1231,10 @@ contains
! apply the second set of givens rotations to a. ! apply the second set of givens rotations to a.
! !
do j = 1, nm1 do j = 1, nm1
if (abs(w(j)) > one) cs = one/w(j) if (abs(w(j)) > 1) cs = 1/w(j)
if (abs(w(j)) > one) sn = sqrt(one-cs**2) if (abs(w(j)) > 1) sn = sqrt(1-cs**2)
if (abs(w(j)) <= one) sn = w(j) if (abs(w(j)) <= 1) sn = w(j)
if (abs(w(j)) <= one) cs = sqrt(one-sn**2) if (abs(w(j)) <= 1) cs = sqrt(1-sn**2)
do i = 1, m do i = 1, m
temp = cs*a(i,j) + sn*a(i,n) temp = cs*a(i,j) + sn*a(i,n)
a(i,n) = -sn*a(i,j) + cs*a(i,n) a(i,n) = -sn*a(i,j) + cs*a(i,n)
@ -1248,7 +1244,7 @@ contains
end subroutine r1mpyq end subroutine r1mpyq
pure subroutine r1updt(m,n,s,ls,u,v,w,sing) pure subroutine r1updt(m,n,s,ls,u,v,w,sing)
use const_and_precisions, only : zero, one, giant=>comp_huge use const_and_precisions, only : giant=>comp_huge
! arguments ! arguments
integer, intent(in) :: m, n, ls integer, intent(in) :: m, n, ls
logical, intent(out) :: sing logical, intent(out) :: sing
@ -1346,8 +1342,8 @@ contains
do nmj = 1, nm1 do nmj = 1, nm1
j = n - nmj j = n - nmj
jj = jj - (m - j + 1) jj = jj - (m - j + 1)
w(j) = zero w(j) = 0
if (v(j) /= zero) then if (v(j) /= 0) then
! !
! determine a givens rotation which eliminates the ! determine a givens rotation which eliminates the
! j-th element of v. ! j-th element of v.
@ -1356,8 +1352,8 @@ contains
cotan = v(n)/v(j) cotan = v(n)/v(j)
sn = p5/sqrt(p25+p25*cotan**2) sn = p5/sqrt(p25+p25*cotan**2)
cs = sn*cotan cs = sn*cotan
tau = one tau = 1
if (abs(cs)*giant > one) tau = one/cs if (abs(cs)*giant > 1) tau = 1/cs
else else
tn = v(j)/v(n) tn = v(j)/v(n)
cs = p5/sqrt(p25+p25*tn**2) cs = p5/sqrt(p25+p25*tn**2)
@ -1393,7 +1389,7 @@ contains
! !
sing = .false. sing = .false.
do j = 1, nm1 do j = 1, nm1
if (w(j) /= zero) then if (w(j) /= 0) then
! !
! determine a givens rotation which eliminates the ! determine a givens rotation which eliminates the
! j-th element of the spike. ! j-th element of the spike.
@ -1402,8 +1398,8 @@ contains
cotan = s(jj)/w(j) cotan = s(jj)/w(j)
sn = p5/sqrt(p25+p25*cotan**2) sn = p5/sqrt(p25+p25*cotan**2)
cs = sn*cotan cs = sn*cotan
tau = one tau = 1
if (abs(cs)*giant > one) tau = one/cs if (abs(cs)*giant > 1) tau = 1/cs
else else
tn = w(j)/s(jj) tn = w(j)/s(jj)
cs = p5/sqrt(p25+p25*tn**2) cs = p5/sqrt(p25+p25*tn**2)
@ -1429,7 +1425,7 @@ contains
! !
! test for zero diagonal elements in the output s. ! test for zero diagonal elements in the output s.
! !
if (s(jj) == zero) sing = .true. if (s(jj) == 0) sing = .true.
jj = jj + (m - j + 1) jj = jj + (m - j + 1)
end do end do
! !
@ -1440,7 +1436,7 @@ contains
s(l) = w(i) s(l) = w(i)
l = l + 1 l = l + 1
end do end do
if (s(jj) == zero) sing = .true. if (s(jj) == 0) sing = .true.
! !
end subroutine r1updt end subroutine r1updt

28
src/vendor/numint.f90 vendored
View File

@ -1,6 +1,6 @@
module numint module numint
use const_and_precisions, only : wp_, zero, one use const_and_precisions, only : wp_
implicit none implicit none
@ -16,10 +16,10 @@ contains
integer :: i integer :: i
real(wp_) :: s0,s1,s2 real(wp_) :: s0,s1,s2
s = zero s = 0
s0 = zero s0 = 0
s1 = zero s1 = 0
s2 = zero s2 = 0
do i = 2, n-1, 2 do i = 2, n-1, 2
s1 = s1+fi(i-1) s1 = s1+fi(i-1)
s0 = s0+fi(i) s0 = s0+fi(i)
@ -39,7 +39,7 @@ contains
real(wp_), intent(out) :: s real(wp_), intent(out) :: s
integer :: i integer :: i
s = zero s = 0
do i = 1, n-1 do i = 1, n-1
s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i)) s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i))
end do end do
@ -110,11 +110,11 @@ contains
! !
! initialize running sums to zero. ! initialize running sums to zero.
! !
flag = zero flag = 0
result = zero result = 0
cor11 = zero cor11 = 0
errest = zero errest = 0
area = zero area = 0
nofun = 0 nofun = 0
if (a .eq. b) return if (a .eq. b) return
! !
@ -124,7 +124,7 @@ contains
nim = 1 nim = 1
x0 = a x0 = a
x(16) = b x(16) = b
qprev = zero qprev = 0
f0 = fun(x0) f0 = fun(x0)
stone = (b - a) / 16.0_wp_ stone = (b - a) / 16.0_wp_
x(8) = (x0 + x(16)) / 2.0_wp_ x(8) = (x0 + x(16)) / 2.0_wp_
@ -174,7 +174,7 @@ contains
! !
! current level is levmax. ! current level is levmax.
! !
flag = flag + one flag = flag + 1
exit exit
end if end if
if (nofun .gt. nofin) then if (nofun .gt. nofin) then
@ -243,7 +243,7 @@ contains
! !
! make sure errest not less than roundoff level. ! make sure errest not less than roundoff level.
! !
if (errest .eq. zero) return if (errest .eq. 0) return
do do
temp = abs(result) + errest temp = abs(result) + errest
if (temp .ne. abs(result)) return if (temp .ne. abs(result)) return