convert_regcoilpm_to_focus

#!/usr/bin/env python

print
print "Usage: " + __file__ + " <regcoil_out.*.nc> <ilambda>"
print "  ilambda is the 0-based index in the lambda scan you want to use."

import sys
import numpy as np

if len(sys.argv) != 3:
    print "Error! You must supply 2 arguments"
    exit(1)

filename = sys.argv[1]
if filename[:12] != 'regcoil_out.':
    print "Error! First argument must be a regcoil_out file"
    exit(1)

try:
    ilambda = int(sys.argv[2])
except:
    print "Error! Unable to convert second argument to an integer."
    exit(1)


print "Opening file "+filename
from scipy.io import netcdf
f = netcdf.netcdf_file(filename,'r',mmap=False)
# We use 'lambdas' instead of 'lambda' to avoid conflict with python's keyword lambda.
nlambda = f.variables['nlambda'][()]
lambdas = f.variables['lambda'][()]
if ilambda > nlambda-1:
    print "Error! Selected ilambda ("+str(ilambda)+") exceeds the number of lambda indices in the regcoil_out file ("+str(nlambda)+") - 1."
    exit(1)

try:
    normal_coil = f.variables['normal_coil'][()]
except:
    print "Error! The field normal_coil was not present in the regcoil_out file. This probably means save_level was set to 3 (the default) in regcoil_pm. Change it to 1."
    exit(1)

chi2_M = f.variables['chi2_M'][()]
chi2_B = f.variables['chi2_B'][()]
print "Using 0-based index ilambda =",ilambda
print "This ilambda corresponds to lambda = {:.3E},  chi^2_B = {:.3E},  chi^2_M = {:.3E}".format(lambdas[ilambda],chi2_B[ilambda],chi2_M[ilambda])

nfp = f.variables['nfp'][()]
sign_normal = f.variables['sign_normal'][()]
d_initial = f.variables['d_initial'][()]
ntheta_coil = f.variables['ntheta_coil'][()]
nzeta_coil = f.variables['nzeta_coil'][()]
ns_magnetization = f.variables['ns_magnetization'][()]
ns_integration = f.variables['ns_integration'][()]
theta_coil = f.variables['theta_coil'][()]
zeta_coil = f.variables['zeta_coil'][()]
zetal_coil = f.variables['zetal_coil'][()]
r_coil = f.variables['r_coil'][()]
s_magnetization = f.variables['s_magnetization'][()]
s_integration = f.variables['s_integration'][()]
s_weights = f.variables['s_weights'][()]
magnetization_vector = f.variables['magnetization_vector'][()]
norm_normal_coil = f.variables['norm_normal_coil'][()]
Jacobian_coil = f.variables['Jacobian_coil'][()]
interpolate_magnetization_to_integration = f.variables['interpolate_magnetization_to_integration'][()]

f.close()

dtheta = theta_coil[1] - theta_coil[0]
dzeta  =  zeta_coil[1] -  zeta_coil[0]

#print "r_coil.shape:",r_coil.shape
#print "interpolate_magnetization_to_integration.shape.shape:",interpolate_magnetization_to_integration.shape
#print "magnetization_vector.shape:",magnetization_vector.shape
#print "Jacobian_coil.shape:",Jacobian_coil.shape

focus_filename = 'focus.regcoil_pm'
f = open(focus_filename,'w')
f.write("# Total number of coils \n")
f.write('{:d} \n'.format(nzeta_coil*ntheta_coil*ns_integration*nfp))
print(ns_integration, ntheta_coil,  nzeta_coil, nfp)

icoil = 0
for js in range(ns_integration):
    for izeta in range(nzeta_coil):
        for itheta in range(ntheta_coil):
            for il in range(nfp):
                izetal = il * nzeta_coil + izeta
                coszeta = np.cos(zetal_coil[izetal])
                sinzeta = np.sin(zetal_coil[izetal])
                icoil += 1
                # Evaluate the position vector in the magnetization region:
                x = r_coil[izetal,itheta,0] + sign_normal * s_integration[js] * d_initial * normal_coil[izetal,itheta,0] / norm_normal_coil[izeta,itheta]
                y = r_coil[izetal,itheta,1] + sign_normal * s_integration[js] * d_initial * normal_coil[izetal,itheta,1] / norm_normal_coil[izeta,itheta]
                z = r_coil[izetal,itheta,2] + sign_normal * s_integration[js] * d_initial * normal_coil[izetal,itheta,2] / norm_normal_coil[izeta,itheta]
                MR = 0
                Mzeta = 0
                MZ = 0
                # Interpolate from the s_magnetization radial grid to the s_integration radial grid:
                for ks in range(ns_magnetization):
                    MR    += magnetization_vector[ilambda,0,ks,izeta,itheta] * interpolate_magnetization_to_integration[ks,js]
                    Mzeta += magnetization_vector[ilambda,1,ks,izeta,itheta] * interpolate_magnetization_to_integration[ks,js]
                    MZ    += magnetization_vector[ilambda,2,ks,izeta,itheta] * interpolate_magnetization_to_integration[ks,js]
                # Convert from (R,zeta) to (x,y) using zeta, not using atan2(y,x), which coincide only at s=0!
                MX = MR * coszeta + Mzeta * (-sinzeta)
                MY = MR * sinzeta + Mzeta * coszeta
                # Convert from magnetization to a discrete dipole
                m_eff_x = dtheta * dzeta * s_weights[js] * Jacobian_coil[js,izeta,izeta] * MX
                m_eff_y = dtheta * dzeta * s_weights[js] * Jacobian_coil[js,izeta,izeta] * MY
                m_eff_z = dtheta * dzeta * s_weights[js] * Jacobian_coil[js,izeta,izeta] * MZ
                moment = np.sqrt(m_eff_x*m_eff_x + m_eff_y*m_eff_y + m_eff_z*m_eff_z)
                # Compute the dipole direction needed by FOCUS
                m_x = m_eff_x / moment
                m_y = m_eff_y / moment
                m_z = m_eff_z / moment
                mt = np.arccos(m_z)
                mp = np.arctan2(m_y, m_x)

                f.write("#-----------------{:10d}---------------------------\n".format(icoil))
                f.write("#coil_type     coil_name\n")
                f.write("  2  pm_{:010d}\n".format(icoil))
                f.write("#  Lc  ox   oy   oz  Ic  I  mt  mp\n")
                f.write("0  {:.15E}  {:.15E}  {:.15E}  0  {:.15E}  {:.15E}  {:.15E}\n".format(x,y,z,moment,mt,mp))

f.close()
print "Conversion complete."