Skip to content
Snippets Groups Projects
Commit 78d68a08 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

implementation of global miller

parent 257a31cc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/geometry_mod.F90
+ 7
− 4
View file @ 78d68a08
...@@ -18,6 +18,8 @@ implicit none ...@@ -18,6 +18,8 @@ implicit none
REAL(xp), PUBLIC, PROTECTED :: s_delta = 0._xp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr REAL(xp), PUBLIC, PROTECTED :: s_delta = 0._xp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr
REAL(xp), PUBLIC, PROTECTED :: zeta = 0._xp ! squareness REAL(xp), PUBLIC, PROTECTED :: zeta = 0._xp ! squareness
REAL(xp), PUBLIC, PROTECTED :: s_zeta = 0._xp ! '' szeta = r dzeta/dr REAL(xp), PUBLIC, PROTECTED :: s_zeta = 0._xp ! '' szeta = r dzeta/dr
REAL(xp), PUBLIC, PROTECTED :: theta1 = 0._xp ! for Miller global
REAL(xp), PUBLIC, PROTECTED :: theta2 = 0._xp !
! to apply shift in the parallel z-BC if shearless ! to apply shift in the parallel z-BC if shearless
REAL(xp), PUBLIC, PROTECTED :: shift_y = 0._xp ! for Arno <3 REAL(xp), PUBLIC, PROTECTED :: shift_y = 0._xp ! for Arno <3
REAL(xp), PUBLIC, PROTECTED :: Npol = 1._xp ! number of poloidal turns (real for 3D zpinch studies) REAL(xp), PUBLIC, PROTECTED :: Npol = 1._xp ! number of poloidal turns (real for 3D zpinch studies)
...@@ -80,7 +82,8 @@ CONTAINS ...@@ -80,7 +82,8 @@ CONTAINS
! Read the input parameters ! Read the input parameters
IMPLICIT NONE IMPLICIT NONE
NAMELIST /GEOMETRY/ geom, q0, shear, eps,& NAMELIST /GEOMETRY/ geom, q0, shear, eps,&
kappa, s_kappa,delta, s_delta, zeta, s_zeta,& ! For miller kappa, s_kappa,delta, s_delta, zeta, s_zeta,& ! For Miller
theta1, theta2,& ! for Miller global
parallel_bc, shift_y, Npol, PB_PHASE parallel_bc, shift_y, Npol, PB_PHASE
READ(lu_in,geometry) READ(lu_in,geometry)
PB_PHASE = .false. PB_PHASE = .false.
...@@ -139,15 +142,15 @@ CONTAINS ...@@ -139,15 +142,15 @@ CONTAINS
kappa = 1._xp kappa = 1._xp
C_y = 1._xp C_y = 1._xp
Cyq0_x0 = 1._xp Cyq0_x0 = 1._xp
CASE('miller','Miller') CASE('miller','Miller','Miller_global','miller_global')
CALL speak('Miller geometry') CALL speak('Miller geometry')
IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for Miller geometry" IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for Miller geometry"
IF(MODULO(FLOOR(Npol),2) .EQ. 0) THEN IF(MODULO(FLOOR(Npol),2) .EQ. 0) THEN
write(*,*) "Npol must be odd for Miller, (Npol = ",Npol,")" write(*,*) "Npol must be odd for Miller, (Npol = ",Npol,")"
ERROR STOP ERROR STOP
ENDIF ENDIF
call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta) call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta,theta1,theta2)
call get_miller(eps,major_R,major_Z,q0,shear,FLOOR(Npol),alpha_MHD,edge_opt,& call get_miller(geom,eps,major_R,major_Z,q0,shear,FLOOR(Npol),alpha_MHD,edge_opt,&
C_y,C_xy,Cyq0_x0,xpdx_pm_geom,gxx,gxy,gxz,gyy,gyz,gzz,& C_y,C_xy,Cyq0_x0,xpdx_pm_geom,gxx,gxy,gxz,gyy,gyz,gzz,&
dBdx,dBdy,dBdz,hatB,jacobian,hatR,hatZ,dxdR,dxdZ) dBdx,dBdy,dBdz,hatB,jacobian,hatR,hatZ,dxdR,dxdZ)
CASE('circular','circ') CASE('circular','circ')
......
src/miller_mod.F90
+ 168
− 78
View file @ 78d68a08
...@@ -21,11 +21,12 @@ MODULE miller ...@@ -21,11 +21,12 @@ MODULE miller
real(xp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta real(xp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
real(xp) :: thetaShift real(xp) :: thetaShift
real(xp) :: thetak, thetad real(xp) :: thetak, thetad
real(xp) :: asurf, theta1, theta2, theta3, delta2, delta3, Raxis, Zaxis
CONTAINS CONTAINS
!>Set defaults for miller parameters !>Set defaults for miller parameters
subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_) subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_,theta1_,theta2_)
real(xp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_ real(xp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_,theta1_,theta2_
rho = -1._xp rho = -1._xp
kappa = kappa_ kappa = kappa_
s_kappa = s_kappa_ s_kappa = s_kappa_
...@@ -36,15 +37,24 @@ CONTAINS ...@@ -36,15 +37,24 @@ CONTAINS
drR = 0._xp drR = 0._xp
drZ = 0._xp drZ = 0._xp
thetak = 0._xp thetak = 0._xp
thetad = 0._xp thetad = 0._xp
asurf = 0.54_xp
theta1 = theta1_
theta2 = theta2_
theta3 = 0._xp
delta2 = 1._xp
delta3 = 1._xp
Raxis = 1._xp
Zaxis = 0._xp
end subroutine set_miller_parameters end subroutine set_miller_parameters
!>Get Miller metric, magnetic field, jacobian etc. !>Get Miller metric, magnetic field, jacobian etc.
subroutine get_miller(trpeps,major_R,major_Z,q0,shat,Npol,amhd,edge_opt,& subroutine get_miller(geom,trpeps,major_R,major_Z,q0,shat,Npol,amhd,edge_opt,&
C_y,C_xy,Cyq0_x0,xpdx_pm_geom,gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,dBdx_,dBdy_,dBdz_,& C_y,C_xy,Cyq0_x0,xpdx_pm_geom,gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,dBdx_,dBdy_,dBdz_,&
Bfield_,jacobian_,R_hat_,Z_hat_,dxdR_,dxdZ_) Bfield_,jacobian_,R_hat_,Z_hat_,dxdR_,dxdZ_)
!!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CHARACTER(len=16), INTENT(IN) :: geom
real(xp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio) real(xp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio)
real(xp), INTENT(INOUT) :: major_R ! major radius real(xp), INTENT(INOUT) :: major_R ! major radius
real(xp), INTENT(INOUT) :: major_Z ! major Z real(xp), INTENT(INOUT) :: major_Z ! major Z
...@@ -73,8 +83,8 @@ CONTAINS ...@@ -73,8 +83,8 @@ CONTAINS
real(xp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong real(xp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong
real(xp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt real(xp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt
real(xp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt real(xp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt
! real(xp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng real(xp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng
! real(xp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng real(xp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng
real(xp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1 real(xp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1
real(xp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr real(xp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr
...@@ -91,9 +101,9 @@ CONTAINS ...@@ -91,9 +101,9 @@ CONTAINS
real(xp), dimension(1:total_nz+ngz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out real(xp), dimension(1:total_nz+ngz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
real(xp), dimension(1:total_nz+ngz):: R_out, Z_out, dxdR_out, dxdZ_out real(xp), dimension(1:total_nz+ngz):: R_out, Z_out, dxdR_out, dxdZ_out
real(xp):: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full real(xp):: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full
!real(xp) :: Lnorm, Psi0 ! currently module-wide defined anyway real(xp) :: Lnorm, Psi0 ! currently module-wide defined anyway
real(xp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz real(xp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz
! real(xp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter real(xp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter
logical:: bMaxShift logical:: bMaxShift
integer:: i, k, iBmax integer:: i, k, iBmax
...@@ -124,76 +134,156 @@ CONTAINS ...@@ -124,76 +134,156 @@ CONTAINS
do while (bMaxShift) do while (bMaxShift)
!flux surface parametrization !flux surface parametrization
thAdj = theta + thetaShift thAdj = theta + thetaShift
if (zeta/=0._xp .or. s_zeta/=0._xp) then SELECT CASE (geom)
R = R0 + rho*Cos(thAdj + d_inv*Sin(thAdj)) CASE('Miller','miller')
Z = Z0 + kappa*rho*Sin(thAdj + zeta*Sin(2*thAdj)) if (zeta/=0._xp .or. s_zeta/=0._xp) then
R = R0 + rho*Cos(thAdj + d_inv*Sin(thAdj))
R_rho = drR + Cos(thAdj + d_inv*Sin(thAdj)) - s_delta*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj)) Z = Z0 + kappa*rho*Sin(thAdj + zeta*Sin(2*thAdj))
Z_rho = drZ + kappa*s_zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
+ kappa*Sin(thAdj + zeta*Sin(2*thAdj)) + kappa*s_kappa*Sin(thAdj + zeta*Sin(2*thAdj)) R_rho = drR + Cos(thAdj + d_inv*Sin(thAdj)) - s_delta*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
Z_rho = drZ + kappa*s_zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
R_theta = -(rho*(1 + d_inv*Cos(thAdj))*Sin(thAdj + d_inv*Sin(thAdj))) + kappa*Sin(thAdj + zeta*Sin(2*thAdj)) + kappa*s_kappa*Sin(thAdj + zeta*Sin(2*thAdj))
Z_theta = kappa*rho*(1 + 2*zeta*Cos(2*thAdj))*Cos(thAdj + zeta*Sin(2*thAdj))
R_theta = -(rho*(1 + d_inv*Cos(thAdj))*Sin(thAdj + d_inv*Sin(thAdj)))
R_theta_theta = -(rho*(1 + d_inv*Cos(thAdj))**2*Cos(thAdj + d_inv*Sin(thAdj))) & Z_theta = kappa*rho*(1 + 2*zeta*Cos(2*thAdj))*Cos(thAdj + zeta*Sin(2*thAdj))
+ d_inv*rho*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
Z_theta_theta = -4*kappa*rho*zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) & R_theta_theta = -(rho*(1 + d_inv*Cos(thAdj))**2*Cos(thAdj + d_inv*Sin(thAdj))) &
- kappa*rho*(1 + 2*zeta*Cos(2*thAdj))**2*Sin(thAdj + zeta*Sin(2*thAdj)) + d_inv*rho*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
else Z_theta_theta = -4*kappa*rho*zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
Rcirc = rho*Cos(thAdj - thetad + thetak) - kappa*rho*(1 + 2*zeta*Cos(2*thAdj))**2*Sin(thAdj + zeta*Sin(2*thAdj))
Zcirc = rho*Sin(thAdj - thetad + thetak) else
Relong = Rcirc Rcirc = rho*Cos(thAdj - thetad + thetak)
Zelong = Zcirc + (-1 + kappa)*rho*Sin(thAdj - thetad + thetak) Zcirc = rho*Sin(thAdj - thetad + thetak)
RelongTilt = Relong*Cos(thetad - thetak) - Zelong*Sin(thetad - thetak) Relong = Rcirc
ZelongTilt = Zelong*Cos(thetad - thetak) + Relong*Sin(thetad - thetak) Zelong = Zcirc + (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
Rtri = RelongTilt - rho*Cos(thAdj) + rho*Cos(thAdj + delta*Sin(thAdj)) RelongTilt = Relong*Cos(thetad - thetak) - Zelong*Sin(thetad - thetak)
Ztri = ZelongTilt ZelongTilt = Zelong*Cos(thetad - thetak) + Relong*Sin(thetad - thetak)
RtriTilt = Rtri*Cos(thetad) + Ztri*Sin(thetad) Rtri = RelongTilt - rho*Cos(thAdj) + rho*Cos(thAdj + delta*Sin(thAdj))
ZtriTilt = Ztri*Cos(thetad) - Rtri*Sin(thetad) Ztri = ZelongTilt
R = R0 + RtriTilt RtriTilt = Rtri*Cos(thetad) + Ztri*Sin(thetad)
Z = Z0 + ZtriTilt ZtriTilt = Ztri*Cos(thetad) - Rtri*Sin(thetad)
R = R0 + RtriTilt
drRcirc = Cos(thAdj - thetad + thetak) Z = Z0 + ZtriTilt
drZcirc = Sin(thAdj - thetad + thetak)
drRelong = drRcirc drRcirc = Cos(thAdj - thetad + thetak)
drZelong = drZcirc - (1 - kappa - kappa*s_kappa)*Sin(thAdj - thetad + thetak) drZcirc = Sin(thAdj - thetad + thetak)
drRelongTilt = drRelong*Cos(thetad - thetak) - drZelong*Sin(thetad - thetak) drRelong = drRcirc
drZelongTilt = drZelong*Cos(thetad - thetak) + drRelong*Sin(thetad - thetak) drZelong = drZcirc - (1 - kappa - kappa*s_kappa)*Sin(thAdj - thetad + thetak)
drRtri = drRelongTilt - Cos(thAdj) + Cos(thAdj + delta*Sin(thAdj)) & drRelongTilt = drRelong*Cos(thetad - thetak) - drZelong*Sin(thetad - thetak)
- s_delta*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj)) drZelongTilt = drZelong*Cos(thetad - thetak) + drRelong*Sin(thetad - thetak)
drZtri = drZelongTilt drRtri = drRelongTilt - Cos(thAdj) + Cos(thAdj + delta*Sin(thAdj)) &
drRtriTilt = drRtri*Cos(thetad) + drZtri*Sin(thetad) - s_delta*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
drZtriTilt = drZtri*Cos(thetad) - drRtri*Sin(thetad) drZtri = drZelongTilt
R_rho = drR + drRtriTilt drRtriTilt = drRtri*Cos(thetad) + drZtri*Sin(thetad)
Z_rho = drZ + drZtriTilt drZtriTilt = drZtri*Cos(thetad) - drRtri*Sin(thetad)
R_rho = drR + drRtriTilt
dtRcirc = -(rho*Sin(thAdj - thetad + thetak)) Z_rho = drZ + drZtriTilt
dtZcirc = rho*Cos(thAdj - thetad + thetak)
dtRelong = dtRcirc dtRcirc = -(rho*Sin(thAdj - thetad + thetak))
dtZelong = dtZcirc + (-1 + kappa)*rho*Cos(thAdj - thetad + thetak) dtZcirc = rho*Cos(thAdj - thetad + thetak)
dtRelongTilt = dtRelong*Cos(thetad - thetak) - dtZelong*Sin(thetad - thetak) dtRelong = dtRcirc
dtZelongTilt = dtZelong*Cos(thetad - thetak) + dtRelong*Sin(thetad - thetak) dtZelong = dtZcirc + (-1 + kappa)*rho*Cos(thAdj - thetad + thetak)
dtRtri = dtRelongTilt + rho*Sin(thAdj) - rho*(1 + delta*Cos(thAdj))*Sin(thAdj + delta*Sin(thAdj)) dtRelongTilt = dtRelong*Cos(thetad - thetak) - dtZelong*Sin(thetad - thetak)
dtZtri = dtZelongTilt dtZelongTilt = dtZelong*Cos(thetad - thetak) + dtRelong*Sin(thetad - thetak)
dtRtriTilt = dtRtri*Cos(thetad) + dtZtri*Sin(thetad) dtRtri = dtRelongTilt + rho*Sin(thAdj) - rho*(1 + delta*Cos(thAdj))*Sin(thAdj + delta*Sin(thAdj))
dtZtriTilt = dtZtri*Cos(thetad) - dtRtri*Sin(thetad) dtZtri = dtZelongTilt
R_theta = dtRtriTilt dtRtriTilt = dtRtri*Cos(thetad) + dtZtri*Sin(thetad)
Z_theta = dtZtriTilt dtZtriTilt = dtZtri*Cos(thetad) - dtRtri*Sin(thetad)
R_theta = dtRtriTilt
dtdtRcirc = -(rho*Cos(thAdj - thetad + thetak)) Z_theta = dtZtriTilt
dtdtZcirc = -(rho*Sin(thAdj - thetad + thetak))
dtdtRelong = dtdtRcirc dtdtRcirc = -(rho*Cos(thAdj - thetad + thetak))
dtdtZelong = dtdtZcirc - (-1 + kappa)*rho*Sin(thAdj - thetad + thetak) dtdtZcirc = -(rho*Sin(thAdj - thetad + thetak))
dtdtRelongTilt = dtdtRelong*Cos(thetad - thetak) - dtdtZelong*Sin(thetad - thetak) dtdtRelong = dtdtRcirc
dtdtZelongTilt = dtdtZelong*Cos(thetad - thetak) + dtdtRelong*Sin(thetad - thetak) dtdtZelong = dtdtZcirc - (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
dtdtRtri = dtdtRelongTilt + rho*Cos(thAdj) - rho*(1 + delta*Cos(thAdj))**2*Cos(thAdj + delta*Sin(thAdj)) & dtdtRelongTilt = dtdtRelong*Cos(thetad - thetak) - dtdtZelong*Sin(thetad - thetak)
+ delta*rho*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj)) dtdtZelongTilt = dtdtZelong*Cos(thetad - thetak) + dtdtRelong*Sin(thetad - thetak)
dtdtZtri = dtdtZelongTilt dtdtRtri = dtdtRelongTilt + rho*Cos(thAdj) - rho*(1 + delta*Cos(thAdj))**2*Cos(thAdj + delta*Sin(thAdj)) &
dtdtRtriTilt = dtdtRtri*Cos(thetad) + dtdtZtri*Sin(thetad) + delta*rho*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
dtdtZtriTilt = dtdtZtri*Cos(thetad) - dtdtRtri*Sin(thetad) dtdtZtri = dtdtZelongTilt
R_theta_theta = dtdtRtriTilt dtdtRtriTilt = dtdtRtri*Cos(thetad) + dtdtZtri*Sin(thetad)
Z_theta_theta = dtdtZtriTilt dtdtZtriTilt = dtdtZtri*Cos(thetad) - dtdtRtri*Sin(thetad)
endif R_theta_theta = dtdtRtriTilt
Z_theta_theta = dtdtZtriTilt
endif
CASE('Miller_global','miller_global')
rho_a = rho/aSurf
CN2 = (-1._xp + Delta2**2)/(1._xp + Delta2**2)
CN3 = (-1._xp + Delta3**2)/(1._xp + Delta3**3)
psiN = rho_a**2 + CN2*rho_a**2 + CN3*rho_a**3
drpsiN = 2._xp*rho_a + 2._xp*CN2*rho_a + 3._xp*CN3*rho_a**2
xAng = 2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 - 27._xp*CN3**2*Cos(3._xp*(thAdj + theta3))**2*psiN
drxAng = -27._xp*CN3**2*Cos(3._xp*(thAdj + theta3))**2*drpsiN
dtxAng = -12._xp*CN2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2*Sin(2._xp*(thAdj + theta2)) &
+ 162._xp*CN3**2*Cos(3._xp*(thAdj + theta3))*psiN*Sin(3._xp*(thAdj + theta3))
dtdtxAng = 486._xp*CN3**2*Cos(6._xp*(thAdj + theta3))*psiN + 24._xp*CN2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2))) &
*(-(Cos(2._xp*(thAdj + theta2))*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))) + 2._xp*CN2*Sin(2._xp*(thAdj + theta2))**2)
yAng = 3._xp*Sqrt(3._xp)*CN3*Cos(3._xp*(thAdj + theta3))*Sqrt(psiN)*Sqrt(2._xp*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 + xAng)
dryAng = (yAng*(drpsiN/psiN + drxAng/(2._xp*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 + xAng)))/2._xp
dtyAng = yAng*(-3._xp*Tan(3._xp*(thAdj + theta3)) &
+ (-12._xp*CN2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2*Sin(2._xp*(thAdj + theta2)) + dtxAng) &
/(2._xp*(2*(1._xp + CN2*Cos(2*(thAdj + theta2)))**3 + xAng)))
dtdtyAng = dtyAng**2/yAng + yAng*(-9._xp/Cos(3._xp*(thAdj + theta3))**2 &
- (-12._xp*CN2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2*Sin(2._xp*(thAdj + theta2)) + dtxAng)**2 &
/(2._xp*(2._xp*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 + xAng)**2) &
+ (-24._xp*CN2*Cos(2._xp*(thAdj + theta2))*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2 &
+ 48._xp*CN2**2*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))*Sin(2._xp*(thAdj + theta2))**2 + dtdtxAng) &
/(2._xp*(2._xp*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**3 + xAng)))
rAng = atan2(yAng,xAng)/3.
drrAng = (-(yAng*drxAng) + xAng*dryAng)/(3._xp*(xAng**2 + yAng**2))
dtrAng = (-(yAng*dtxAng) + xAng*dtyAng)/(3._xp*(xAng**2 + yAng**2))
dtdtrAng = (-6._xp*dtrAng**2*yAng)/xAng &
+ ((2._xp*yAng*dtxAng**2)/xAng - 2._xp*dtxAng*dtyAng - yAng*dtdtxAng + xAng*dtdtyAng)/(3._xp*(xAng**2 + yAng**2))
rShape = (aSurf*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))/Cos(3._xp*(thAdj + theta3))*&
&(-1._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng)))/(3._xp*CN3)
drrShape = (aSurf*(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))/Cos(3._xp*(thAdj + theta3))*&
&(Sqrt(3._xp)*Cos(rAng) - Sin(rAng))*drrAng)/(3._xp*CN3)
dtrShape = rShape*((-2._xp*CN2*Sin(2._xp*(thAdj + theta2)))/(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))&
&+ 3._xp*Tan(3._xp*(thAdj + theta3)) &
&+ ((Sqrt(3._xp)*Cos(rAng) - Sin(rAng))*dtrAng)/(-1._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng)))
dtdtrShape = dtrShape**2/rShape + rShape*(9._xp - (4._xp*CN2*Cos(2._xp*(thAdj + theta2)))/(1._xp + CN2*Cos(2._xp*(thAdj + theta2))) &
- (4._xp*CN2**2*Sin (2._xp*(thAdj + theta2))**2)/(1._xp + CN2*Cos(2._xp*(thAdj + theta2)))**2 + 9._xp*Tan (3._xp*(thAdj + theta3))**2 &
+ ((-4._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng))*dtrAng**2)/(-1._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng))**2 &
+ ((Sqrt(3._xp)*Cos(rAng) - Sin(rAng))*dtdtrAng)/(-1._xp + Cos(rAng) + Sqrt(3._xp)*Sin(rAng)))
do i=1._xp,np
if (abs(CN3*cos(3._xp*(thAdj(i)+theta3)))<1e-8) then
rShape(i) = aSurf*Sqrt(psiN/(1._xp + CN2*Cos(2*(thAdj(i) + theta2))))
drrShape(i) = (rShape(i)*drpsiN)/(2._xp*psiN)
dtrShape(i) = (CN2*rShape(i)*Sin(2._xp*(thAdj(i) + theta2)))/(1._xp + CN2*Cos(2._xp*(thAdj(i) + theta2)))
dtdtrShape(i) = dtrShape(i)**2/rShape(i) &
+ (2*(CN2**2 + CN2*Cos(2._xp*(thAdj(i) + theta2)))*rShape(i))/(1._xp + CN2*Cos(2._xp*(thAdj(i) + theta2)))**2
endif
if (abs(1._xp + CN2*Cos(2._xp*(thAdj(i) + theta2)))<1e-8) then
rShape(i) = aSurf*Sqrt(psiN)
drrShape(i) = 1._xp
dtrShape(i) = 0._xp
dtdtrShape(i) = 0._xp
endif
enddo
Rcenter = Raxis - (-R0 + Raxis)*rho_a**2
Zcenter = Zaxis - rho_a**2*(-Z0 + Zaxis)
drRcenter = -2._xp*(-R0 + Raxis)*rho_a
drZcenter = -2._xp*rho_a*(-Z0 + Zaxis)
R = Rcenter + Cos(thAdj)*rShape
Z = rShape*Sin(thAdj) + Zcenter
R_rho = (drRcenter + Cos(thAdj)*drrShape)/aSurf ! Adjust for rho deriv, not rho_a deriv
Z_rho = (drZcenter + Sin(thAdj)*drrShape)/aSurf ! Adjust for rho deriv, not rho_a deriv
R_theta = -(rShape*Sin(thAdj)) + Cos(thAdj)*dtrShape
Z_theta = Cos(thAdj)*rShape + Sin(thAdj)*dtrShape
R_theta_theta = -(Cos(thAdj)*rShape) - 2._xp*Sin(thAdj)*dtrShape + Cos(thAdj)*dtdtrShape
Z_theta_theta = -(rShape*Sin(thAdj)) + 2._xp*Cos(thAdj)*dtrShape + Sin(thAdj)*dtdtrShape
END SELECT
!dl/dtheta !dl/dtheta
dlp=(R_theta**2+Z_theta**2)**0.5_xp dlp=(R_theta**2+Z_theta**2)**0.5_xp
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment