Vmec conversion

[RalfMackenbach]

Hello,

I added VMEC conversion scripts, right now it's standalone. It includes a test case at the end of the script. Bref is currently based on psi_edge, so not the same as B0. Definition of alpha may also be different. Let me know if I can help

Cheers

[gsnoep]

As discussed yesterday, the quickest way to integrate this PR is to write a reader for the Equilibrium class that puts the VMEC information directly into Equilibrium.fluxsurfaces, e.g.:

def read_vmec(self,f_path=None,fit_modes=2**8):
        
    # read the g-file
    with io.netcdf_file(f_path,'r',mmap=False) as file:
        vmec = file.variables
    
    # retrieve rmnc, zmns, xn, ns, and nfp
    rmnc = vmec['rmnc']
    zmns = vmec['zmns']
    xn = (vmec['xn']).data
    xm = (vmec['xm']).data
    ns = (vmec['ns']).data
    nfp = (vmec['nfp']).data
    # retrieve pres and iotas
    pres = vmec['presf'].data
    iota = vmec['iotaf'].data
    pres = pres[1::]
    iota = iota[1::]
    phi_last = vmec['phi'].data[-1]
    edge_toroidal_flux_over_2pi = phi_last / (2*np.pi)
    L_reference = vmec['Aminor_p'].data

    # make grid of theta and zeta
    nzeta= 1
    theta = np.linspace(0,2*np.pi,num=fit_modes)
    zeta = np.linspace(0,2*np.pi/nfp,num=nzeta,endpoint=False)
    # number of modes
    nmodes = len(xn)
    iradius = ns-1
    R = np.zeros((ns-1,fit_modes,nzeta))
    Z = np.zeros((ns-1,fit_modes,nzeta))
    for iradius in range(ns-1):
        for itheta in range(fit_modes):
            for izeta in range(nzeta):
                for imode in range(nmodes):
                    angle = xm[imode]*theta[itheta] - xn[imode]*zeta[izeta]
                    R[iradius,itheta,izeta] = R[iradius,itheta,izeta] + rmnc[iradius,imode]*math.cos(angle)
                    Z[iradius,itheta,izeta] = Z[iradius,itheta,izeta] + zmns[iradius,imode]*math.sin(angle)
                
    self.fluxsurfaces['R'] = R
    self.fluxsurfaces['Z'] = Z

    # add other radial profiles for phi, psi, rho_tor, rho_pol, q, B_tor, etc here
        
    return

and then we can adjust the LocalEquilibrium class to allow for supplying an Equilibrium that already contains flux surface information to fit.