calculate_vmec_geometry Subroutine

public subroutine calculate_vmec_geometry(vmec_filename, nalpha, alpha0, nzgrid, zeta_center, rectangular_cross_section, number_of_field_periods_inputfile, s_inputfile, vmec_surface_option, verbose, s, safety_factor_q, shat, L_reference, B_reference, nfp_out, sign_toroidal_flux, alpha, zeta, bmag, b_dot_grad_zeta, grad_alpha_grad_alpha, grad_alpha_grad_psit, grad_psit_grad_psit, gds23_psitalpha, gds24_psitalpha, gds25_psitalpha, gds26_psitalpha, gbdrift_alpha, gbdrift0_psit, cvdrift_alpha, cvdrift0_psit, theta, B_sub_zeta, B_sub_theta, psit_displacement_fac, gradzeta_gradpsit_R2overB2, gradzeta_gradalpha_R2overB2, b_dot_grad_zeta_RR, ierr)

             CALCULATE TRIPLE PRODUCTS OF B, ∇s, ∇α                !

   CALCULATE THE CARTESIAN COMPONENTS OF THE GRADIENT VECTORS       !

Use R(ialpha,izeta) and Z(ialpha,izeta), to compute X = R * cos(zeta) Y = R * sin(zeta)

We can write the magnetic field in covariant representation (the components transform like the transformation of the reference axes) B = B_s ∇s + B_theta ∇θ + B_zeta ∇ζ Where be assume that B_s = 0 for a vacuum magnetic field B = B_theta ∇θ + B_zeta ∇ζ

We can also write the magnetic field in contravariant representation (the components transform inversely to the transformation of the reference axes) B = (1/sqrt(g)) (dpsi_t/ds) (iota * e_theta + e_zeta) So we have a tangent basis (∇s, ∇θ, ∇ζ) and the dual basis (e_s, e_theta, e_zeta) ∇s . e_s = 1 ∇θ . e_theta = 1 ∇ζ . e_zeta = 1 ∇s . e_theta = 0 ∇s . e_zeta = 0 ∇θ . e_s = 0 ... e_s = (dX/ds) e_X + (dY/ds) e_Y + (dZ/ds) e_Z e_ζ = (dX/dζ) e_X + (dY/dζ) e_Y + (dZ/dζ) e_Z e_θ = (dX/dθ) e_X + (dY/dθ) e_Y + (dZ/dθ) e_Z

The two bases are related through the dual relations sqrt(g) = 1 / (∇s . ∇θ x ∇ζ)
e_i = sqrt(g) (e^j x e^k) e_s = sqrt(g) (∇θ x ∇ζ) e_ζ = sqrt(g) (∇s x ∇θ) e_θ = sqrt(g) (∇ζ x ∇s) ∇s = 1/sqrt(g) (e_θ x e_ζ) ∇ζ = 1/sqrt(g) (e_s x e_θ) ∇θ = 1/sqrt(g) (e_ζ x e_s)

If we multiply the covariant representation with the contravariant representation we get B^2 = (1/sqrt(g)) (dpsi_t/ds) (iotaB_theta + B_zeta) The Jacobian sqrt(g) of the magnetic coordinate system can thus be calculated as sqrt(g) = (dpsi_t/ds) (iotaB_theta + B_zeta)/B^2

Use the dual relations to get the Cartesian components of ∇s, ∇θ, and ∇ζ ∇s = 1/sqrt(g) (e_θ x e_ζ) = 1/sqrt(g) ((dX/dθ)e_X + (dY/dθ)e_Y + (dZ/dθ)e_Z) x ((dX/dζ)e_X + (dY/dζ)e_Y + (dZ/dζ)e_Z) ∇s . e_X = 1/sqrt(g) ( dY/dθ * dZ/dζ - dZ/dθ * dY/dζ ) ∇s . e_Y = 1/sqrt(g) ( dZ/dθ * dX/dζ - dX/dθ * dZ/dζ ) ∇s . e_Z = 1/sqrt(g) ( dX/dθ * dY/dζ - dY/dθ * dX/dζ )

        PERFORM THE SINE/COSINE FOURIER TRANSFORMATIONS             !

We constructed the arrays and along the chosen field lines, so now we can perform the sine/consine fourier transformations along these field lines.

Note that R, Z, and Lambda use the (xm,xn) mode numbers, while all the other quantities use the (xm_nyq,xn_nyq) mode numbers. B and Lambda are on the half mesh, so their radial derivatives are on the full mesh. R and Z are on the full mesh, so their radial derivatives are on the half mesh.

For the derivatives with respect to and we have, e.g., R = sum_{imn} R_{imn}cos(angle) dR/dtheta = sum_{imn} d(R_{imn}cos(angle))/dtheta = - m * sum_{imn} R_{imn}sin(angle) dR/dzeta = sum_{imn} d(R_{imn}cos(angle))/dzeta = n * sum_{imn} R_{imn}*sin(angle)

                 QUANTITIES ON THE FLUX SURFACE                     !

Evaluate several radial-profile functions at the flux surface we ended up choosing in the routine. Here , , and and is a module variables

                    Choose the flux surface                         !

Determine which flux surface to use, based on and . Possible values of vmec_surface_option: 0 = Use the exact radius requested. 1 = Use the nearest value of the VMEC half grid. 2 = Use the nearest value of the VMEC full grid. Here is a module variable, and , , , , are subroutine variables.

              Print information of the output file                  !

                    Test the input variables                        !

                      Test the VMEC arrays                          !

Do some sanity checking to ensure the VMEC arrays have some expected properties. Here , , , , , are module variables

                          Test iota                              !

If the conversion to has been done correctly, we should find that iota = (B dot grad theta_pest) / (B dot grad zeta)

                           Test ∇ζ                               !

We have the following relations X = Rcos(ζ) ζ = arccos(X/R) Y = Rsin(ζ) ζ = arcsin(Y/R) R^2 = X^2 + Y^2 Therefore, we should have dζ/dX = d(arccos(X/R))/dX = -Y/R^2 = -sin(ζ)/R
dζ/dY = d(arcsin(Y/R))/dX = X/R^2 = cos(ζ)/R
dζ/dZ = 0

                        Test Jacobian                            !

sqrt(g) = 1 / (∇s . ∇θ x ∇ζ) = dX/ds dY/dθ dZ/dζ + dY/ds dZ/dθ dX/dζ + dZ/ds dX/dθ dY/dζ - dX/ds dZ/dθ dY/dζ - dY/ds dX/dθ dZ/dζ - dZ/ds dY/dθ dX/dζ 1/sqrt(g) = (∇s . ∇θ x ∇ζ) = (∇s.e_X) (∇θ.e_Y) (∇ζ.e_Z) + (∇s.e_Y) (∇θ.e_Z) (∇ζ.e_X) + (∇s.e_Z) (∇θ.e_X) (∇ζ.e_Y) - (∇s.e_X) (∇θ.e_Z) (∇ζ.e_Y) + (∇s.e_Y) (∇θ.e_X) (∇ζ.e_Z) + (∇s.e_Z) (∇θ.e_Y) (∇ζ.e_X)

       Test covariant and contravariant componets of B           !

We have the following dual basis e_s = (dX/ds) e_X + (dY/ds) e_Y + (dZ/ds) e_Z e_ζ = (dX/dζ) e_X + (dY/dζ) e_Y + (dZ/dζ) e_Z e_θ = (dX/dθ) e_X + (dY/dθ) e_Y + (dZ/dθ) e_Z Therefore, the covariant components of the magnetic field are B_sub_s = B . e_s = B_X * dX/ds + B_Y * dY/ds + B_Z * dZ/ds B_sub_zeta = B . e_ζ = B_X * dX/dzeta + B_Y * dY/dzeta + B_Z * dZ/dzeta B_sub_theta = B . e_θ = B_X * dX/dtheta + B_Y * dY/dtheta + B_Z * dZ/dtheta And the contravariant components are B_sup_s = B . ∇s = (B.e_X)(∇s.e_X) + (B.e_Y)(∇s.e_Y) + (B.e_Z)(∇s.e_Z) B_sup_zeta = B . ∇ζ = (B.e_X)(∇ζ.e_X) + (B.e_Y)(∇ζ.e_Y) + (B.e_Z)(∇ζ.e_Z) B_sup_theta = B . ∇θ = (B.e_X)(∇θ.e_X) + (B.e_Y)(∇θ.e_Y) + (B.e_Z)(∇θ.e_Z)

                  Check that two arrays match                       !

Verify that array1 = array2 to within a relative tolerance. Verify that |array1| = 0 to within a relative tolerance.

                   Initialize geometry arrays                       !

                   Deallocate geometry arrays                       !

Arguments

                        Output quantities                           !

                       Internal variables                           !