# =====================================================
#      PAS_Moments.py
#   Read PAS data, read PAS calibration, calculate moments 
#   Versions
#   2020-04-21 V.1.0, Fedorov, initial
#   2020-04-21 V.1.1, Fedorov, Density
#   2020-04-22 V.1.2, Fedorov, Velocity
#   2020-04-26 V.2.0, Fedorov, Package pas_process using
#   2020-04-26 V.2.1, Fedorov, Temperature
#   2020-05-01 V.2.2, Fedorov, All calibtaions in the paskage
#   2020-05-01 V.2.3, Fedorov, odd-even elevation tables
# ====================================================

# ----------------------------------------------------
#  Packages
# ----------------------------------------------------
import math
import numpy
from io import StringIO

# ----------------------------------------------------
# CDF reader
# ---------------------------------------------------
import spacepy
import datetime
import time
from spacepy import pycdf

#------------------------------------------------------
# PAS package
#------------------------------------------------------
import pas_process as pas

#-----------------------------------------------------
# Debugging
#-----------------------------------------------------
# import pdb

# -----------------------------------------------------
# Data Intefaces
# -----------------------------------------------------

dDir = "../DATA_L1/"
# normalFile = "20200415/solo_L1_swa-pas-3dn_20200415_V00.cdf"
# normalFile = "20200424/solo_L1_swa-pas-3dn_20200424_V00.cdf"
# normalFile = "20200428/solo_L1_swa-pas-3dn_20200428_V00.cdf"
normalFile = "20200507/solo_L1_swa-pas-3d_20200507_V00.cdf"
# normalFile = "20200506/solo_L1_swa-pas-3d_20200506_V00.cdf"

#------------------------------------------
# Read data file in the range of the file
# -----------------------------------------
recordN = 0
normalCDF = pycdf.CDF(dDir+normalFile)
cdfSize = normalCDF['Epoch'].shape[0]
# print("CDF Size = ",cdfSize)
firstepoch = normalCDF['Epoch'][0]
lastepoch = datetime.datetime(firstepoch.year,firstepoch.month,firstepoch.day,firstepoch.hour)
isolastepoch = '{:.19}'.format(lastepoch.isoformat())

# for recordN in range(cdfSize):
for recordN in range(1030,1041):
    epoch = normalCDF['Epoch'][recordN]
    isoepoch = '{:.19}'.format(epoch.isoformat())
    # pdb.set_trace() # STOP for debugging

    spectrum3D = normalCDF['COUNTS'][recordN,:,:,:] # [ CEMs, Elevarion, Energies]
    sE = normalCDF['FIRST_ENERGY'][recordN]
    sEl = normalCDF['FIRST_ELEVATION'][recordN]
    nE = normalCDF['NB_ENERGY'][recordN]
    nEl = normalCDF['NB_ELEVATION'][recordN]

    # --------------------------------------------
    # Calculate moments
    #---------------------------------------------
    moments = pas.calcMoments(spectrum3D, sE, sEl, nE, nEl)    

    # ---------------------------------------------
    # density
    # ---------------------------------------------
    N = moments["N"]
    
    # ---------------------------------------------
    # Velocities
    # ---------------------------------------------
    Vx = moments["V"][0]
    Vy = moments["V"][1]
    Vz = moments["V"][2]
        
    # ----------------------------------------------
    # Temperature
# ----------------------------------------------
    T = moments["T"]
    Tx = moments["TxTyTz"][0]
    Ty = moments["TxTyTz"][1]
    Tz = moments["TxTyTz"][2]
    
    # -------------------------------
    # print only if the second is changed
    # -------------------------------
    if (isolastepoch != isoepoch):
        print('{}{:4d}{:4d}{:4d}{:4d}{:8.2f}{:9.2f}{:9.2f}{:9.2f}{:9.2f}{:9.2f}{:9.2f}{:9.2f}'.format(
            isoepoch, sE,nE,sEl,nEl,N, Vx*1.e-5,Vy*1.e-5,Vz*1.e-5,T,Tx,Ty,Tz))
        isolastepoch = isoepoch
# ----------------------------------------------------

#-------------------------------------------
# For DPU, print the clean array
#-------------------------------------------
sampleFile = open("SamplleSpectrum", "w")
#-------------------------------------------
# Legende print
#-------------------------------------------
sampleFile.write(
    '# RealPAScount V.2.0, 11 May 2020\n'+
    '# Vx = {:8.1f} km/s'.format(Vx*1.e-5)+' Vx = {:8.1f} km/s'.format(Vy*1.e-5)+' Vz = {:8.1f} km/s\n'.format(Vz*1.e-5) +
    '# T  = {:8.2f} eV\n'.format(T) + 
    '# N =  {:8.2f} cm-3\n'.format(N) + 
    '# Format:\n'+
    '# Energy step number\n'+
    '# Elev. 0, count(CEM 0), count(CEM 1), ... count(CEM 10)\n'+
    '# Elev. 1, count(CEM 0), count(CEM 1), ... count(CEM 10)\n'+
    '# Elev. ..........................\n'+
    '# Elev. 8, count(CEM 0), count(CEM 1), ... count(CEM 10)\n'+
    '#----------------------------------------------------------\n')
#for ie in range(96):
#    sampleFile.write('{:3d}\n'.format(ie))
#    for iel in range(9):
#        sampleFile.write('{:3d}'.format(iel))
#        for iaz in range(11):
#            sampleFile.write('{:10.0f}'.format(spectrum3Dclean[iaz,iel,ie]))
#        sampleFile.write('\n')
sampleFile.close()