Returns fraction of gyroangles that can produce a reaction with given inputs
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(LocalEMFields), | intent(in) | :: | fields | Electromagneticfields in beam coordinates |
||
| real(kind=Float64), | intent(in) | :: | E3 | E3 charged fusion product energy [keV] |
||
| real(kind=Float64), | intent(in) | :: | E1 | E1 fast-ion energy [keV] |
||
| real(kind=Float64), | intent(in) | :: | pitch | pitch fast ion pitch relative to the field |
||
| type(LocalProfiles), | intent(in) | :: | plasma | Plasma Paramters in beam coordinates |
||
| real(kind=Float64), | intent(in), | dimension(3) | :: | v3_xyz | Charged fusion product velocity in beam coorindates |
|
| real(kind=Float64), | intent(out) | :: | pgyro | pgyro DeltaE_3*\partial\gam/\partial E_3/pi |
||
| real(kind=Float64), | intent(out) | :: | gam0 | Gyro angle of fast ion [rad] |
subroutine get_pgyro(fields,E3,E1,pitch,plasma,v3_xyz,pgyro,gam0)
!+ Returns fraction of gyroangles that can produce a reaction with
!+ given inputs
type(LocalEMFields), intent(in) :: fields
!+ Electromagneticfields in beam coordinates
real(Float64), intent(in) :: E3
!+ E3 charged fusion product energy [keV]
real(Float64), intent(in) :: E1
!+ E1 fast-ion energy [keV]
real(Float64), intent(in) :: pitch
!+ pitch fast ion pitch relative to the field
type(LocalProfiles), intent(in) :: plasma
!+ Plasma Paramters in beam coordinates
real(Float64), dimension(3), intent(in) :: v3_xyz
!+ Charged fusion product velocity in beam coorindates
real(Float64), intent(out) :: pgyro
!+ pgyro DeltaE_3*\partial\gam/\partial E_3/pi
real(Float64), intent(out) :: gam0
!+ Gyro angle of fast ion [rad]
real(Float64), dimension(3) :: a_hat, vrot
real(Float64) :: JMeV,JkeV,mp,Q,norm_v3,norm_v1,vpar,vperp,vb,va,E3max,E3min
real(Float64) :: phip,cosphip,cosphib,cosphia,rhs,gammaplus,gammaminus
real(Float64) :: eps,bracket,ccminus,ccplus,bbminus,bbplus,v3minus,v3plus
real(Float64) :: E3minus,E3plus,v3plusmag,v3minusmag,v3maxmag,v3minmag
real(Float64) :: DeltaE3
pgyro = 0.d0
gam0 = 0.d0
! Preliminaries [SI units]
JMeV = 1.60218d-13 ! Conversion factor from MeV to Joules
JkeV = 1.60218d-16 ! Conversion factor from keV to Joules
mp = H1_amu*mass_u ![kg]
Q = 4.04*JMeV ![J]
norm_v3 = sqrt(E3/(H1_amu*v2_to_E_per_amu)) / 100 ![m/s]
norm_v1 = sqrt(E1/(beam_mass*v2_to_E_per_amu)) / 100 ![m/s]
vpar = norm_v1*pitch ![m/s]
vperp = norm_v1*sqrt(1-pitch**2) ![m/s]
vrot = plasma%vrot/100.d0 ![m/s]
a_hat = cross_product(fields%b_norm, v3_xyz) / norm2(v3_xyz) !(b,a,c)
!! First, check E3 limits are valid for the given inputs
!! Assumes cos(gamma) = +/- 1 and vrot = 0
cosphip = dot_product(v3_xyz, fields%b_norm) / norm2(v3_xyz)
phip = acos(cosphip)
!! Get plasma rotation components in (b,a,c) coords
if (all(vrot.eq.0.d0)) then
vb = 0.d0
va = 0.d0
else
cosphib = dot_product(vrot, fields%b_norm) / norm2(vrot)
vb = norm2(vrot)*cosphib ![m/s]
cosphia = dot_product(vrot, a_hat) / norm2(vrot)
va = norm2(vrot)*cosphia ![m/s]
endif
! LHS coefficient teeny (step #0)
bracket = vperp*(sin(phip)-2*va/norm_v3)
if (abs(bracket).lt.1.d-5) then !'Case 0'
return
endif
! Find E3 limits for these parameters (step #1)
ccminus = 1.5*Q/mp + 0.5*norm_v1**2 - 2*vpar*vb + 0.5*norm2(vrot) - vperp*va
ccplus = 1.5*Q/mp + 0.5*norm_v1**2 - 2*vpar*vb + 0.5*norm2(vrot) + vperp*va
bbminus = -vperp*sin(phip) - (vpar+vb)*cos(phip) - va*sin(phip)
bbplus = vperp*sin(phip) - (vpar+vb)*cos(phip) - va*sin(phip)
v3minus = 0.5*(-bbminus + sqrt(bbminus**2+4*ccminus))
v3plus = 0.5*(-bbplus + sqrt(bbplus**2+4*ccplus))
E3minus = 0.5*mp*v3minus**2 / JkeV
E3plus = 0.5*mp*v3plus**2 / JkeV
E3min = min(E3plus,E3minus) ![keV]
E3max = max(E3plus,E3minus) ![keV]
! Now E3plus and E3minus are charged fusion product energies at edge of bin
! 'Case 1'
E3plus = E3 + 0.5*ctable%dE
E3minus = E3 - 0.5*ctable%dE
! Bin misses gamma curve altogether (step #2)
if ((E3plus.le.E3min) .or. (E3minus.ge.E3max)) then !'Case 2'
return
endif
! Get gammaplus and gammaminus (step #3)
v3plusmag = sqrt(2*E3plus*JkeV/mp)
v3minusmag = sqrt(2*E3minus*JkeV/mp)
v3maxmag = sqrt(2*E3plus*JkeV/mp)
v3minmag = sqrt(2*E3minus*JkeV/mp)
if (E3plus.gt.E3max) then !'Case 3'
norm_v3 = v3maxmag
rhs = norm_v3 - 1.5*Q/(mp*norm_v3) - (vpar+vb)*cos(phip) - va*sin(phip) &
- 0.5*(norm_v1**2 + norm2(vrot))/norm_v3 + 2*vpar*vb/norm_v3
if (abs(rhs).gt.abs(bracket)) then
if (rhs*bracket.gt.0.) then
gammaplus = 0.
else
gammaplus = pi
endif
else
gammaplus = acos(rhs/bracket)
endif
else
norm_v3 = v3plusmag
rhs = norm_v3 - 1.5*Q/(mp*norm_v3) - (vpar+vb)*cos(phip) - va*sin(phip) &
- 0.5*(norm_v1**2 + norm2(vrot))/norm_v3 + 2*vpar*vb/norm_v3
if (abs(rhs).gt.abs(bracket)) then
if (rhs*bracket.gt.0.) then
gammaplus = 0.
else
gammaplus = pi
endif
else
gammaplus = acos(rhs/bracket)
endif
endif
if (E3minus.lt.E3min) then !'Case 3'
norm_v3 = v3minmag
rhs = norm_v3 - 1.5*Q/(mp*norm_v3) - (vpar+vb)*cos(phip) - va*sin(phip) &
- 0.5*(norm_v1**2 + norm2(vrot))/norm_v3 + 2*vpar*vb/norm_v3
if (abs(rhs).gt.abs(bracket)) then
if (rhs*bracket.gt.0.) then
gammaminus = 0.
else
gammaminus = pi
endif
else
gammaminus = acos(rhs/bracket)
endif
else
norm_v3 = v3minusmag
rhs = norm_v3 - 1.5*Q/(mp*norm_v3) - (vpar+vb)*cos(phip) - va*sin(phip) &
- 0.5*(norm_v1**2 + norm2(vrot))/norm_v3 + 2*vpar*vb/norm_v3
if (abs(rhs).gt.abs(bracket)) then
if (rhs*bracket.gt.0.) then
gammaminus = 0.
else
gammaminus = pi
endif
else
gammaminus = acos(rhs/bracket)
endif
endif
gam0 = (gammaplus + gammaminus) / 2.d0 ! return mean gyroangle
!!! should there be modifications in gamma in if statements below
if ((E3plus.ge.E3max) .and. (E3minus.le.E3min)) then !'Case 4a'
pgyro = 1.d0
return
else if ((E3plus.lt.E3max) .and. (E3minus.gt.E3min)) then !'Case 4b'
pgyro = abs(gammaplus-gammaminus)/pi
return
else if((E3plus.ge.E3max) .and. (E3minus.gt.E3min)) then
!'Case 4c'
if (gammaplus.gt.gammaminus) then
pgyro = (pi-gammaminus)/pi
return
else
pgyro = gammaminus/pi
return
endif
else if ((E3plus.lt.E3max) .and. (E3minus.le.E3min)) then
!'Case 4d'
if (gammaminus.gt.gammaplus) then
pgyro = (pi-gammaplus)/pi
return
else
pgyro = gammaplus/pi
return
endif
endif
endsubroutine get_pgyro