apps120811_flux_recalc.regenerate_fluxfile

''' Regenerate the precipitation flux

Different algorithm can provide different precipitation flux.
There are indeed several ways, which I examined.

:mod:`precipitation`
    Precipitation map considering the thermal velocity.
    Integral over pitch angle >90.  The plasma parameter
    at 1.2 Rm is mapped along the B field to the surface.
    A data file is produced as "MHD_PrecipitationFlux_YF_footprint.dat".

:mod:`precipitation_net`
    Precipitation map using the net flux, i.e. N x V, along B.
    This is the one Xianzhe calculated first.
    A data file is NOT produces, since it will give the same
    (or similar) results of "MHD_PrecipitationFlux_Jia.dat".

:mod:`precipitation_notrace`
    Precipitation map considering the thermal velocity.
    Integral over pitch angle >90.  The plasma parameter
    at 1.2 Rm is mapped to the nadir.
    A data file is produced as "MHD_PrecipitationFlux_YF_nadir.dat".

If those produced files are placed to the path specified in
the path specified in ``[pep]`` section's ``mhddatapath`` entry,
they can be read from the :mod:`pyana.pep.mhddata` module.
'''

from pyana.pep import mhddata

from . import precipitation
from . import precipitation_notrace


def mk_footprint(mf, pp):
    retval = precipitation.main(mf=mf, pp=pp, dlat=1., dlon=1.)

    lons = retval[2]
    lats = retval[3]
    flx = retval[4]   # Flux shape of (lon, lat)

    fp = open('MHD_PrecipitationFlux_YF_footprint.dat', 'w')
    print('VARIABLES="longitude(deg)", "latitude(deg)", "ParticleFlux(/cm2/s)"', file=fp)
    for ilon in range(len(lons)):
        for ilat in range(len(lats)):
            print(lons[ilon], lats[ilat], flx[ilon, ilat], file=fp)

    fp.close()
            
def mk_nadir(mf, pp):
    retval = precipitation_notrace.main(mf=mf, pp=pp, dlat=1., dlon=1.)

    lons = retval[2]
    lats = retval[3]
    flx = retval[4]   # Flux shape of (lon, lat)

    fp = open('MHD_PrecipitationFlux_YF_nadir.dat', 'w')
    print('VARIABLES="longitude(deg)", "latitude(deg)", "ParticleFlux(/cm2/s)"', file=fp)
    for ilon in range(len(lons)):
        for ilat in range(len(lats)):
            print(lons[ilon], lats[ilat], flx[ilon, ilat], file=fp)

    fp.close()
            

def main():
    mf = mhddata.MagneticField1201()
    mf.set_interpolate_strategy('project')
    pp = mhddata.PlasmaParameter1205()

    mk_footprint(mf, pp)
    mk_nadir(mf, pp)

if __name__ == "__main__":
    main()