# ================================================================= # pas_process.py # PAS procession procedures: # - calibration files read # - spectra corrections # - disiribution function # Versions: # 26 Apr 2020 . V.1.0, Fedorov, The first package Version # 30 Apr 2020 . V.1.1, Fedorov, Single count remove # 01 May 2020 . V.1.2, Fedorov, All calirations files read # 13 May 2020 . V.1.3, Fedorov, Simple energy table also # 15 May 2020 . V.1.4, Fedorov, elevation correction table read, GF table read # 15 May 2020 . V.1.5, Fedorov, All tables read during import, all tables are properties of the package # 16 May 2020. V.1.6, Fedorov, Moments calculation are here # 17 May 2020. V.1.7, Fedorov, VDF # 17 May 2020. V.1.8, Fedorov, correct pressure tensot inJ/cm3 # ================================================================== import numpy import math from io import StringIO # -------------------------------------- # For debugging # ------------------------------------- # import pdb # pdb.set_trace() # ----------------------------------------- # Global variables corresponding for internal use # ----------------------------------------- calFilesDir = "../../CALIBRATION/DOCLOG/PFM_CALIBRATION_RECORD/" catalogFile = "solo_CAL_catalog_20200514_V05.txt" catalogDict = { # catalog Dictionaty "EVS": "", # Simple energy table "ELF": "", # Full Elavation Table "AZF": "", # Full Azimuth table "GFF": "", # GF table (* accumulation time ) "CNF": "", # Conversion to dN ( * accumulation time ) "YIL": "", # Ghost correction "ELD": "" # Elevation Odd-even perturbance } eTable = {"emin":[], "emax":[], "ecent":[], "vcent":[]} elTable = {"elmin":[], "elmax":[], "elcent":[], "elsin": [], "elcos":[], "index":[]} azTable = {"azmin":[], "azmax":[], "azcent":[], "azsin":[], "azcos":[], "index":[]} elCorrTab = {"deg":[], "rad":[]} yeildArr = numpy.zeros((11,9,11,9),dtype = float) # [ Az_BIN, El_BIN, Az_Cell, El_Cell] gfTable = numpy.zeros((11,9),dtype = float) # GF table (* accumulation time ) cNArr = numpy.zeros((11,9,96),dtype = float) # CN vXarr = numpy.zeros((11,9,96), dtype = float) # Vx vYarr = numpy.zeros((11,9,96), dtype = float) # Vy vZarr = numpy.zeros((11,9,96), dtype = float) # Vz moments = {"N":0.0, "V":[0.0,0.0,0.0], "PxxPyyPzzPxyPxzPyz":[0.0,0.0,0.0,0.0,0.0,0], "T":0.0, "TxTyTz":[0.0,0.0,0.0]} vdf = numpy.zeros((11,9,96),dtype = float) # 3D velocity distribution function P2T = 5.2197e-13 # recalculation from Pressure (in mass = 1) to T P2Jcm3 = 1.6726e-31 # conversion the pressure for internal use J/cm3 #------------------------------------------ # Function to read calibration catalog # ----------------------------------------- def catalogRead(): cFile = open(calFilesDir+catalogFile, "r") fullList = cFile.readlines() listSize = len(fullList) for it in range(13,listSize): litems = fullList[it].split() # ------------------------------- # fined the index in the dictionary # ------------------------------- if litems[1] in catalogDict: catalogDict[litems[1]] = litems[2] cFile.close() # ------------------------------------------------- #------------------------------------------ # Function to read calibration SIMPLE Energy table # eArr is centers # ----------------------------------------- def eTableRead(): eArrMin = [] eArrMax = [] eFile = open(calFilesDir+catalogDict["EVS"], "r") for il in range(8): line = eFile.readline() # legend lines # print(line) for il in range(96): line = eFile.readline() # data lines lener = line.split() eArrMin.append(float(lener[1])) eArrMax.append(float(lener[2])) eFile.close() eMinArr = numpy.array(eArrMin,dtype = float) eMaxArr = numpy.array(eArrMax,dtype = float) eArr = numpy.sqrt(eMinArr * eMaxArr) vArr = 13.841e5*numpy.sqrt(eArr) #cm/s #--------------------------------------------------- # For DPU. Print V table #--------------------------------------------------- # for ie in range(96): # print('{:4d}{:8.2f}'.format(ie,vArr[ie]*1.0e-5)) eTable["emin"] = eMinArr eTable["emax"] = eMaxArr eTable["ecent"] = eArr eTable["vcent"] = vArr # return(eTable) # ----------------------------------------- #------------------------------------------ # Function to read calibration FULL Elevation table # elArr is centers # ----------------------------------------- def elTableRead(): elArrMin = numpy.zeros((11,9),dtype = float) elArrMax = numpy.zeros((11,9),dtype = float) elArr = numpy.zeros((11,9),dtype = float) elsArr = numpy.zeros((11,9),dtype = float) # sin elcArr = numpy.zeros((11,9),dtype = float) # cos elFile = open(calFilesDir+catalogDict["ELF"], "r") for il in range(9): line = elFile.readline() # legend lines for il in range(11): # Azimuth line = elFile.readline() # Az number azbin = line.split() iaz = int(azbin[2]) for iel in range(9): line = elFile.readline() lelev = line.split() elArrMin[iaz,iel] = float(lelev[1]) elArrMax[iaz,iel] = float(lelev[2]) if ((elArrMin[iaz,iel] > -400.0) and (elArrMax[iaz,iel] > -400.0)): elArr[iaz,iel] = (elArrMin[iaz,iel] + elArrMax[iaz,iel])/2.0 else : elArr[iaz,iel] = -1.00E+31 elsArr[iaz,iel] = float(lelev[3]) elcArr[iaz,iel] = float(lelev[4]) elFile.close() # print(elArr) # ------------------------------------------ # index of correct values in the table # ------------------------------------------ elCalIndx = numpy.nonzero([elArr > -400.0]) elTable["elmin"] = elArrMin elTable["elmax"] = elArrMax elTable["elcent"] = elArr elTable["elsin"] = elsArr elTable["elcos"] = elcArr elTable["index"] = elCalIndx # print(elArr[elCalIndx]) # elTable = {"elmin": elArrMin, "elmax":elArrMax, "elcent":elArr, # "elsin": elsArr, "elcos":elcArr, "index":elCalIndx} # return(elTable) # ----------------------------------------- #------------------------------------------ # Function to read calibration FULL Azimuth table # azArr is centers # ----------------------------------------- def azTableRead(): azArrMin = numpy.zeros((11,9),dtype = float) azArrMax = numpy.zeros((11,9),dtype = float) azsArr = numpy.zeros((11,9),dtype = float) # sin azcArr = numpy.zeros((11,9),dtype = float) # cos azArr = numpy.zeros((11,9),dtype = float) azFile = open(calFilesDir+catalogDict["AZF"], "r") for il in range(9): line = azFile.readline() # legend lines for ilel in range(9): line = azFile.readline() # elevation number elbin = line.split() iel = int(elbin[2]) for iaz in range(11): line = azFile.readline() # data line lelev = line.split() azArrMin[iaz,iel] = float(lelev[1]) azArrMax[iaz,iel] = float(lelev[2]) if ((azArrMin[iaz,iel] > -400.0) and (azArrMax[iaz,iel] > -400.0)): azArr[iaz,iel] = (azArrMin[iaz,iel] + azArrMax[iaz,iel])/2.0 else : azArr[iaz,iel] = -1.00E+31 azsArr[iaz,iel] = float(lelev[3]) azcArr[iaz,iel] = float(lelev[4]) azFile.close() azCalIndx = numpy.nonzero([azArr > -400.0]) azTable["azmin"] = azArrMin azTable["azmax"] = azArrMax azTable["azcent"] = azArr azTable["azsin"] = azsArr azTable["azcos"] = azcArr azTable["index"] = azCalIndx # azTable = {"azmin": azArrMin, "azmax":azArrMax, "azcent":azArr, # "azsin": azsArr, "azcos":azcArr, "index":azCalIndx} # return(azTable) # ----------------------------------------- #------------------------------------------ # Function to read calibration CN table # ----------------------------------------- def cnTableRead(): index = numpy.arange(1, 12, 1) cnFile = open(calFilesDir+catalogDict["CNF"], "r") for il in range(9): line = cnFile.readline() # legend lines for ie in range(96): line = cnFile.readline() ie = int(line) for iel in range(9): line = cnFile.readline() buff = numpy.genfromtxt(StringIO(line), delimiter=(5,10,10,10,10,10,10,10,10,10,10,10)) ncline = buff[index] indx = numpy.nonzero([ncline >= 0.0]) cNArr[indx[1],iel,ie] = ncline[indx[1]] cnFile.close() # return(cNArr) # --------------------------------------------- #------------------------------------------ # Function to read calibration GF table # ----------------------------------------- def gfTableRead(): index = numpy.arange(1, 12, 1) gfFile = open(calFilesDir+catalogDict["GFF"], "r") for il in range(8): line = gfFile.readline() # legend lines for iel in range(9): line = gfFile.readline() # elevation lines buff = numpy.genfromtxt(StringIO(line), delimiter=(5,11,11,11,11,11,11,11,11,11,11,11)) gfline = buff[index] gfTable[:,iel] = gfline gfFile.close() # --------------------------------------------- #------------------------------------------ # Function to read elevation correction # ----------------------------------------- def elCorrRead(): elCorrArrD = numpy.zeros((96),dtype = float) # Elevation Correction Degres elCorrArrR = numpy.zeros((96),dtype = float) # Elevation Correction Rad elCorrFile = open(calFilesDir+catalogDict["ELD"], "r") for il in range(8): line = elCorrFile.readline() # legend lines for ie in range(96): line = elCorrFile.readline() # calibration line lines = line.split() elCorrArrD[ie] = float(lines[1]) elCorrArrR[ie] = float(lines[2]) elCorrFile.close() elCorrTab["deg"] = elCorrArrD elCorrTab["rad"] = elCorrArrR # --------------------------------------------- # ---------------------------------------------------- # Function to read the ghost yield array # ----------------------------------------------------- def yeildRead(): # it is the string for the file name yeildFile = open(calFilesDir+catalogDict["YIL"], "r") for il in range(10): line = yeildFile.readline() # legend lines for ibazl in range(11): line = yeildFile.readline() # BIN azimuthal number lines = line.split() ibaz = int(lines[2]) for ibell in range(9): line = yeildFile.readline() # BIN elevation number lines = line.split() ibel = int(lines[2]) for icazl in range(11): line = yeildFile.readline() # line of one CAL azimuth lines = line.split() icaz = int(lines[2]) yeildArr[ibaz,ibel,icaz,:] = numpy.array(lines[3:12],dtype = float) yeildFile.close() # ----------------------------------------------------------- # ------------------------------------------------------------ # Function correcting one spectrum # ----------------------------------------------------------- def specCorrection(spectrum3D, sEl, nEl, sE, nE): spectrum3DC = numpy.zeros((11,nEl,nE),dtype = float) for ibaz in range(11): for ibel in range(nEl): for ie in range(nE): # --------------------------- # remove single INSULATED count # the summ oh the nex bins shall be greather than # --------------------------- if (spectrum3D[ibaz,ibel,ie] == 1.0): around = 0.0 if (ibaz > 0): around += spectrum3D[ibaz-1,ibel,ie] if (ibaz < 10): around += spectrum3D[ibaz+1,ibel,ie] if (ibel > 0): around += spectrum3D[ibaz,ibel-1,ie] if (ibel < nEl-1): around += spectrum3D[ibaz,ibel+1,ie] if (ie > 0): around += spectrum3D[ibaz,ibel,ie-1] if (ie < nE-1): around += spectrum3D[ibaz,ibel,ie+1] if (around < 2.0): spectrum3D[ibaz,ibel,ie] = 0.0 # ------------------------------------- # remove ghost counts # ------------------------------------ partRemove = spectrum3D[ibaz, ibel, ie] - numpy.sum(spectrum3D[:,:,ie]*yeildArr[ibaz, ibel+sEl,:,sEl:sEl+nEl]) if (partRemove < 0.0): partRemove = 0.0 # print(ibaz,ibel,ie, partRemove,) spectrum3DC[ibaz, ibel, ie] = partRemove return(spectrum3DC) # ----------------------------------------------------------- # ------------------------------------------------------------ # Functions calculating velocities #------------------------------------------------------------- def getVelocities(): elsArr = elTable["elsin"] # sin elcArr = elTable["elcos"] # cos azsArr = azTable["azsin"] # sin azcArr = azTable["azcos"] elCorrR = elCorrTab["rad"] vArr = eTable["vcent"] for ie in range(96): elcArrC = elcArr - elsArr*elCorrR[ie] elsArrC = elsArr + elcArr*elCorrR[ie] vXarr[:,:,ie] = -vArr[ie]*elcArrC*azcArr vYarr[:,:,ie] = vArr[ie]*elcArr*azsArr vZarr[:,:,ie] = -vArr[ie]*elsArrC # ----------------------------------------------------------- # ------------------------------------------------------------ # Functions calculating Velocity Distribution Function #------------------------------------------------------------- def calcVDF(spectrum, sE, sEl, nE, nEl): spectrum3D = numpy.array(spectrum,dtype = float) spectrum3Dnul = numpy.zeros(spectrum3D.shape,dtype = float) vArr4 = numpy.zeros(96,dtype = float) for ie in range(96): vArr4[ie] = math.pow(eTable["vcent"][ie],4) # --------------------------------------------- # Set Nulls out of boundaries # --------------------------------------------- for iel in range(nEl): for ie in range(nE): spectrum3Dnul[:,iel+sEl,ie+sE] = spectrum3D[:,iel+sEl,ie+sE] # ------------------------------------------ # Clean the Spectrum # ------------------------------------------ spectrum3Dclean = specCorrection(spectrum3Dnul, 0, 9, 0, 96) # pdb.set_trace() # ------------------------------------------ # Transform to VDF # ------------------------------------------ for iaz in range(11): for iel in range(sEl,sEl+nEl): if (gfTable[iaz, iel] > 1.0e-12): vdf[iaz,iel,:] = spectrum3Dclean[iaz,iel,:]/(vArr4 * gfTable[iaz, iel]) return(vdf) # ----------------------------------------------------------- # ------------------------------------------------------------ # Functions calculating moments #------------------------------------------------------------- def calcMoments(spectrum, sE, sEl, nE, nEl): spectrum3D = numpy.array(spectrum,dtype = float) spectrum3Dnul = numpy.zeros(spectrum3D.shape,dtype = float) # --------------------------------------------- # Set Nulls out of boundaries # --------------------------------------------- for iel in range(nEl): for ie in range(nE): spectrum3Dnul[:,iel+sEl,ie+sE] = spectrum3D[:,iel+sEl,ie+sE] #------------------------------------------ # Clean the Spectrum #------------------------------------------ spectrum3Dclean = specCorrection(spectrum3Dnul, 0, 9, 0, 96) # --------------------------------------------- # Calculate the density # --------------------------------------------- moments["N"] = 0.0 moments["V"] = [0.0,0.0,0.0] moments["PxxPyyPzzPxyPxzPyz"] = [0.0,0.0,0.0,0.0,0.0,0] moments["T"] = 0.0 moments["TxTyTz"] = [0.0,0.0,0.0] partN = spectrum3Dclean*cNArr N = numpy.sum(partN) moments["N"] = N # --------------------------------------------- # Calculate the velocity # --------------------------------------------- if (N > 0.001): Vx = numpy.sum(partN*vXarr)/N Vy = numpy.sum(partN*vYarr)/N Vz = numpy.sum(partN*vZarr)/N moments["V"] = [Vx,Vy,Vz] # ---------------------------------------------- # Pressure thensor and Temperature # ---------------------------------------------- if (N > 0.001): Pxx = numpy.sum((vXarr[:,:,:] -Vx)*(vXarr[:,:,:] - Vx)*partN) Pyy = numpy.sum((vYarr[:,:,:] -Vy)*(vYarr[:,:,:] - Vy)*partN) Pzz = numpy.sum((vZarr[:,:,:] -Vz)*(vZarr[:,:,:] - Vz)*partN) moments["PxxPyyPzzPxyPxzPyz"] = [Pxx*P2Jcm3, Pyy*P2Jcm3, Pzz*P2Jcm3, numpy.sum((vXarr[:,:,:] -Vx)*(vYarr[:,:,:] - Vy)*partN)*P2Jcm3, numpy.sum((vXarr[:,:,:] -Vx)*(vZarr[:,:,:] - Vz)*partN)*P2Jcm3, numpy.sum((vYarr[:,:,:] -Vy)*(vZarr[:,:,:] - Vz)*partN)*P2Jcm3] moments["T"] = (Pxx + Pyy + Pzz)/N*P2T*2.0/3.0 moments["TxTyTz"] = [Pxx/N*P2T*2.0, Pyy/N*P2T*2.0, Pzz/N*P2T*2.0] return(moments) # ---------------------------------------------------------- # Init Part that reads all calibration files # ---------------------------------------------------------- catalogRead() eTableRead() elTableRead() azTableRead() cnTableRead() elCorrRead() yeildRead() getVelocities() gfTableRead() # ----------------------------------------------------------