https://github.com/SimonUnterstrasser/ColumnModel
Tip revision: 8c4e8c105fdfa552938c3dcbe71de99be72243f3 authored by SimonUnterstrasser on 15 September 2020, 15:26:46 UTC
Update README.md
Update README.md
Tip revision: 8c4e8c1
PlotResults.gcc.py
'''
main plot module that controls the generation of a-posteriori plot from existing simulation data. Simulation data is read (by the use of submodule InputOutput) and plots are generated (by the use of submodule PlotSim) as specified in the caller script
'''
import sys
import os
import math
import numpy as np
import InputOutput as IO
import PlotSim as PS
#GCCif (BIN == 1)
import Referenzloesung as REF
#GCCendif /* (BIN == 1) */
text=None
label=None
title=None
ilastGVonly=None
iDmean=0
iPlotextras = 0
#GCCinclude "params_plot.txt"
#GCCif (((IREF == 3) && (KERNEL != 3)) || ((IREF == 2) && (KERNEL != 2)))
print('check preprocessor settings')
print('wrong combination of IREF and KERNEL')
exit(1)
#GCCendif
#GCCif (((IREF == 3) && (DISCRETE == 0)) || ((IREF == 2) && (DISCRETE == 0)))
print('check preprocessor settings')
print('wrong combination of IREF and DISCRETE')
exit(1)
#GCCendif
#GCCif ((MOM_meanTD == 2) && (COLUMN == 0))
print('check preprocessor settings')
print('wrong combination of MOM_meanTD and COLUMN')
exit(1)
#GCCendif
#GCCif (IMODE == 1)
def plotSingleSimulationGV(fp_in):
#Einlesen der SIP-Daten
nEK_sip_plot, mEK_sip_plot, zEK_sip_plot, nr_SIPs_plot, nr_SIPs_prof, t_vec_GVplot, V, skal_m = IO.readSingleSimulationSIPdata(fp_in)
nz=1
#GCCif (COLUMN == 1)
nz=IO.get_nz(fp_in)
#GCCendif /*(COLUMN == 1)*/
#Erstelle GV-Plot
PS.PlotGV(nEK_sip_plot/nz, mEK_sip_plot, nr_SIPs_plot, t_vec_GVplot, ilastGVonly = ilastGVonly)
#<<<<<<<< end routine plotSingleSimulationGV(fp)
#GCCif (DISCRETE >= 1)
def plotSingleSimulationRStD(fp_in):
#Read SIP data
nEK_sip_plot, mEK_sip_plot,zDummy, nr_SIPs_plot, nr_SIPs_prof, t_vec_GVplot, V, skal_m = IO.readSingleSimulationSIPdata(fp_in)
#Generate plot of relative standard deviation
PS.PlotRelStdDev(mEK_sip_plot, t_vec_GVplot)
#<<<<<<<< end routine plotSingleSimulationRStD(fp)
#GCCendif /* (DISCRETE >= 1) */
def plotSingleSimulationMOM(fp_in,iprof=0,itot=1):
# read Moment data of box model or column model data
# for the box model the design of the plot is straightforward, i.e. temporal evolution of moments
# for the colum model several options exists:
# 1. temporal evolution of column integrated moments (itot = 1)
# 2. profiles of moments at selected times (iprof = 1)
MOMsave, t_vec_MOMsave, skal_m, dV=IO.readSingleSimulationMoments(fp_in)
iSIPplot = 0
nr_SIPs_prof = None
t_vec_GVplot = None
nr_SIPs_mean = None
nr_SIPs_prof_mean = None
#GCCif (addSIPplot > 0)
iSIPplot = 1
nEK_sip_plot, mEK_sip_plot,zDummy, nr_SIPs_plot, nr_SIPs_prof, t_vec_GVplot, V, skal_m = IO.readSingleSimulationSIPdata(fp_in)
nr_SIPs_mean = np.mean(nr_SIPs_plot, axis=0)
nr_SIPs_prof_mean = np.mean(nr_SIPs_prof, axis=0)
#GCCendif /* (addSIPplot > 0) */
#GCCif (STD == 0)
iplotStdDev = 0
#GCCendif /* (STD == 0) */
#GCCif (STD == 1)
iplotStdDev = 1
#GCCendif /* (STD == 1) */
#GCCif (STD == 2)
iplotStdDev = 2
#GCCendif /* (STD == 2) */
nz = 1 # for COLUMN = 0
#GCCif (COLUMN == 1)
#in column model version: 1. call profile plot 2. integrate over column 3. call total moments plot
# >>>> 1. call profile plot
dz=np.zeros(1)
dz[0] = IO.get_dz(fp_in)
nz=np.zeros(1, dtype=int)
nz[0] = IO.get_nz(fp_in)
if (iprof == 1):
if (iSIPplot == 1):
#GCCif (addSIPplot == 1)
nr_SIPs_prof_mean = nr_SIPs_prof_mean*100/dz[0]
#GCCendif /* (addSIPplot == 1) */
#GCCif (addSIPplot == 2)
nr_SIPs_prof_mean = nr_SIPs_prof_mean
#GCCendif /* (addSIPplot == 2) */
pass
for Times_elem in Times_list:
PS.PlotMomentsProf(MOMsave, t_vec_MOMsave, iplot_onlyMOMmean=1,
nz_vec=nz, dz_vec=dz, iplotStdDev=iplotStdDev, iDmean=iDmean,
iTimes_select=Times_elem,
iSIPplot=iSIPplot, nr_SIPs_prof_mean=nr_SIPs_prof_mean, t_vec_GVplot=t_vec_GVplot)
# >>>> 2. integrate over column
MOMsave=np.mean(MOMsave,axis=2)
#GCCendif /*(COLUMN == 1)*/
#GCCif (MOM_meanTD == 2)
MOMsave = MOMsave * dz[0] * nz[0]
#GCCendif /* (MOM_meanTD == 2) */
#GCCif (MOM_meanTD == 3)
MOMsave = MOMsave * dV * nz[0]
#GCCendif /* (MOM_meanTD == 3) */
#Erstelle Momenten-Plot
if (itot == 1):
#GCCif (MOM_skal_m == 0)
PS.PlotMoments(MOMsave,t_vec_MOMsave,iplot_onlyMOMmean=1,iplotStdDev=iplotStdDev,iDmean=iDmean,
iSIPplot=iSIPplot,nr_SIPs_mean=nr_SIPs_mean,t_vec_GVplot=t_vec_GVplot)
#GCCelse
PS.PlotMoments(MOMsave,t_vec_MOMsave,iplot_onlyMOMmean=1,iplotStdDev=iplotStdDev,iDmean=iDmean,
skal_m=skal_m,
iSIPplot=iSIPplot,nr_SIPs_mean=nr_SIPs_mean,t_vec_GVplot=t_vec_GVplot)
#GCCendif /* else (MOM_skal_m == 0) */
#<<<<<<<< end routine plotSingleSimulationMOM(fp)
def plotSingleSimulationFLUX(fp_in):
#GCCif (STD == 0)
iplotStdDev = 0
#GCCendif /* (STD == 0) */
#GCCif (STD == 1)
iplotStdDev = 1
#GCCendif /* (STD == 1) */
#GCCif (STD == 2)
iplotStdDev = 2
#GCCendif /* (STD == 2) */
#read flux data
FluxIn,FluxOut,FluxInAcc,FluxOutAcc,t_vec_GVplot=IO.readSingleSimulationFluxes(fp_in)
#generate flux plot
PS.PlotFluxesTime(FluxIn,FluxOut,FluxInAcc,FluxOutAcc,t_vec_GVplot[1:],iplot_onlyMean=1,iplotStdDev=iplotStdDev)
#<<<<<<<< end routine plotSingleSimulationFLUX(fp_in)
def plotSingleSimulationGVout(fp_in,infalling=0):
# Read SIP data
# returns a list of SIPs of length nr_SIPs_out with nu and m and the time tEK_sip_out each SIP crossed the lower BC
# set infalling = 1 for additional tracking of infalling SIPs
# A is the basal area of the column
nEK_sip_out, mEK_sip_out, tEK_sip_out, nr_SIPs_out, Tsim, A = IO.readSingleSimulationSIP_outfalling_data(fp_in,infalling)
PS.PlotGV(nEK_sip_out/A, mEK_sip_out, nr_SIPs_out, [Tsim], outfallingGV = 1 + infalling)
#<<<<<<<< end routine plotSingleSimulationGV(fp)
#GCCendif /* (IMODE == 1) */
#GCCif (IMODE == 2)
import Misc as FK
import PlotSim as PS
def plotMultipleSimulationsMOM(filepath_input,iTime=0,iProf=0,iFinal=0,iTimesCross=0):
nr_sims=len(filepath_input)
print('call of plotMultipleSimulationMOM')
print('iTime, iProf, iFinal', iTime, iProf, iFinal)
print('nr_sims:',nr_sims)
#if (label is not None) and (len(label) != nr_sims+iaddsimREF):
#print('wrong number of labels are given')
#print('label',label)
#print('nr_sims, iaddsimREF: ', nr_sims, iaddsimREF)
iplot_mode=np.zeros(nr_sims,dtype=int) +2
iplot_mode[0]=1
iplot_mode[-1]=3
nr_Mom = 4
#GCCif (STD == 0)
iplotStdDev = 0
MOM_StdDev = None
MOMprof_StdDev = None
#GCCendif /* (STD == 0) */
#GCCif (STD == 1)
iplotStdDev = 1
#GCCendif /* (STD == 1) */
#GCCif (STD == 2)
iplotStdDev = 2
#GCCendif /* (STD == 2) */
nt_vec = np.zeros(nr_sims,dtype='int')
nr_inst_vec = np.zeros(nr_sims,dtype='int')
dV_vec = np.zeros(nr_sims)
nz_vec = np.zeros(nr_sims,dtype='int')
nt_eval = np.zeros(nr_sims,dtype='int')
#GCCif (COLUMN == 1)
dz_vec =np.zeros(nr_sims)
#find out first the maximum dimensions of the profile data
for i in range(nr_sims):
fp=filepath_input[i]
print('i, fp', i, fp)
#GCCif (BIN != 1)
nz,dz,Tsim,nr_inst,LWC,r0,xf0,ikernel,i_init_1D,i_process=IO.get_MetaData(fp)
nz_vec[i] = nz
dz_vec[i] = dz
nr_inst_vec[i] = nr_inst
#####GCCelse
####ntREF,nzREF,dzREF,dtREF,TsimREF,nbinREF,scalREF,dlnrREF,rgridREF,mgridREF=REF.get_RefMetaData(fp=fp)
####nz_vec[i] = nzREF
####dz_vec[i] = dzREF
####nt_vec[i] = ntREF
####nr_inst_vec[i] = 1
#GCCendif /* (BIN != 1) */
if (iProf == 1):
nz_max = nz_vec.max()
dz_max = dz_vec.max()
nr_inst_max = nr_inst_vec.max()
nt_max=300
MOMsave=np.zeros([nr_sims,nt_max,nz_max,nr_Mom])
t_vec_MOM_vec=np.zeros([nr_sims,nt_max])
if (iplotStdDev == 1):
MOMprof_StdDev = np.zeros([1,nr_sims,nt_max,nz_max,nr_Mom])
if (iplotStdDev == 2):
MOMprof_StdDev = np.zeros([2,nr_sims,nt_max,nz_max,nr_Mom])
if (iFinal == 1):
MOMFinal = np.zeros([nr_sims,nr_Mom])
if (iTimesCross == 1):
TimesCross = np.zeros([nr_sims,2,3]) #second dimension: 2 components for lambda_0 and lambda_2, third dimension for three threshold levels
iSIPplot = 0
nr_SIPs_prof_mean = None
nr_SIPs_mean = None
t_vec_GVplot_vec = None
t_vec_GVplotsingleSim = None
ntGV_vec = None
#GCCif (addSIPplot > 0)
iSIPplot = 1
nz_max = nz_vec.max()
nr_SIPs_mean, nr_SIPs_prof_mean, t_vec_GVplot_vec, ntGV_vec = plotMultipleSimulationGV(filepath_input, igetSIPdataonly=1, nz_max=nz_max, dz_vec=dz_vec) # with igetSIPdataonly=1 the routine does not generate a plot, but returns the SIP data of all simulations
t_vec_GVplotsingleSim = np.squeeze(t_vec_GVplot_vec[0,:])
#GCCendif /* (addSIPplot > 0) */
#GCCendif /*(COLUMN == 1)*/
for i in range (nr_sims):
fp=filepath_input[i]
print('process simulation:',fp)
#GCCif (BIN != 1)
iFinalData = iFinal * 2
if (iFinal == 1) and (iProf == 0) and (iTime == 0):
#faster execution possible as file MomentsMean.dat is processed
MOMmeanFinalsingleSim, dV, skal_m, iFinalData = IO.readSingleSimulationMomentsMeanFinal(fp,iPlotextras=iPlotextras)
print('FinalData: ', i,iFinalData,MOMmeanFinalsingleSim)
MOMFinal[i,:] = MOMmeanFinalsingleSim
# if iFinalData = 1 then reading the MeanMoments data was succesful
# iFinalData =1 can only happen for the above combination of (iFinal = 1) and (iProf = 0) and (iTime = 0)
#iFinalData:
# 0 no iFinal=1 analysis
# 1 iFinal=1 analysis with mean data succesful. Moreover iProf = 0 and iTime = 0 and no further processing of full data is necessary
# 2 iFinal=1 analysis with mean data not succesful. Read normal data.
if (iFinalData != 1):
MOMsingleSim, t_vec_MOMsingleSim, skal_m, dV =IO.readSingleSimulationMoments(fp)
nt_eval[i] = t_vec_MOMsingleSim.size
dV_vec[i] = dV
if (iTimesCross == 1):
#faster execution possible as file MomentsMean.dat is processed
TimesCross_singleSim = IO.readSingleSimulationMomentsMeanFinal(fp,iPlotextras=iPlotextras,iTimesCross=1,thresh=thresholds_TimeCross)
print(i, TimesCross_singleSim)
TimesCross[i,:,:] = TimesCross_singleSim
#GCCelse
ntREF,nzREF,dzREF,dtREF,TsimREF,nbinREF,scalREF,dlnrREF,rgridREF,mgridREF=REF.get_RefMetaData(fp=fp)
MOMsingleSim = REF.get_RefProfileData(nzREF, ntREF, fp, iswap=1, fp=fp)
t_vec_MOMsingleSim = TsimREF/(ntREF-1)*np.arange(ntREF)
nt_eval[i] = t_vec_MOMsingleSim.size
if "Bott" in fp:
nt_eval[i] = ((nt_eval[i]-1)/2)+1
dV_vec[i] = 1.0
skal_m = 1
iFinalData = 2 # ="analysis with mean data not succesful. Read normal data." as for bin data no shortcut exists
print('TsimREF, ntREF:', TsimREF, ntREF)
#GCCendif /* (BIN != 1) */
#GCCif (COLUMN == 0)
nz_vec[i] = 1
#GCCendif /* (COLUMN == 0) */
if (iProf == 1):
nt_vec[i] = t_vec_MOMsingleSim.size
print('AA',nt_vec[i])
if (nt_vec[i] > nt_max):
print ('nt_vec[i], nt_max', nt_vec[i], nt_max)
sys.exit('increase nt_max!!')
#GCCif (STD > 0)
MOMsave[i,:nt_vec[i],:nz_vec[i],:], MOMprof_StdDev[:,i,:nt_vec[i],:nz_vec[i],:] = \
FK.CIO_MOMmean(data_in=MOMsingleSim,ikeep1D=1,igetSTD=iplotStdDev)
#GCCelse
MOMsave[i,:nt_vec[i],:nz_vec[i],:]=FK.CIO_MOMmean(data_in=MOMsingleSim,ikeep1D=1)
#GCCendif /* else (STD > 0) */
t_vec_MOM_vec[i,:nt_vec[i]]=t_vec_MOMsingleSim
if (iTime == 1) or (iFinalData == 2):
#produce plot MeanMoments over time
#individual call for each simulation
iexcludetop = 0
#GCCif (STD > 0)
MOMmean, MOM_StdDev = FK.CIO_MOMmean(data_in=MOMsingleSim,igetSTD=iplotStdDev,iexcludetop=iexcludetop)
#GCCelse
print(MOMsingleSim.shape)
MOMmean=FK.CIO_MOMmean(data_in=MOMsingleSim,iexcludetop=iexcludetop)
print(MOMmean.shape)
#GCCendif /* else (STD > 0) */
if (iFinal == 1):
MOMFinal[i,:]=MOMmean[nt_eval[i]-1,:]
print('AA: ', nt_eval[i]-1,MOMFinal[i,:])
print('sim: ', i)
print('initial moments ', MOMmean[0])
print('skal_m: ', skal_m)
#GCCif (MOM_meanTD == 2)
MOMmean = MOMmean * dz_vec[i] * nz_vec[i]
#GCCif (STD > 0)
MOM_StdDev = MOM_StdDev * dz_vec[i] * nz_vec[i]
#GCCendif /* (STD > 0) */
#GCCendif /* (MOM_meanTD == 2) */
#GCCif (MOM_meanTD == 3)
MOMmean = MOMmean * dV_vec[i] * nz_vec[i]
#GCCif (STD > 0)
MOM_StdDev = MOM_StdDev * dV_vec[i] * nz_vec[i]
#GCCendif /* (STD > 0) */
#GCCendif /* (MOM_meanTD == 3) */
nr_SIPs_mean_pick = None
#GCCif (addSIPplot > 0)
nr_SIPs_mean_pick = nr_SIPs_mean[i,:]
#GCCendif /* (addSIPplot > 0) */
if (iTime == 1):
#GCCif (MOM_skal_m == 0)
PS.PlotMoments(MOMmean,t_vec_MOMsingleSim,iMean=1,iplot_mode=iplot_mode[i],
label=label,title=title,nr_sims=nr_sims,text=text,
iplotStdDev=iplotStdDev, MOM_StdDev=MOM_StdDev, iDmean=iDmean,
iSIPplot=iSIPplot,nr_SIPs_mean=nr_SIPs_mean_pick,t_vec_GVplot=t_vec_GVplotsingleSim,
iexcludetop=iexcludetop)
#GCCelse
PS.PlotMoments(MOMmean,t_vec_MOMsingleSim,iMean=1,iplot_mode=iplot_mode[i],
label=label,title=title,nr_sims=nr_sims,text=text,
iplotStdDev=iplotStdDev, MOM_StdDev=MOM_StdDev, iDmean=iDmean, skal_m=skal_m,
iSIPplot=iSIPplot,nr_SIPs_mean=nr_SIPs_mean_pick,t_vec_GVplot=t_vec_GVplotsingleSim,
iexcludetop=iexcludetop)
#GCCendif /* else (MOM_skal_m == 0) */
#produce plot MomentProfiles
if (iProf == 1):
for Times_elem in Times_list:
PS.PlotMomentsProf(MOMsave, t_vec_MOM_vec, iMultipleSims=1, iMean=1,
nr_inst_vec=nr_inst_vec, nz_vec=nz_vec, nt_vec=nt_vec, dz_vec = dz_vec,
label=label,title=title,text=text,
iplot_onlyMOMmean=1,iTimes_select=Times_elem,
iplotStdDev=iplotStdDev, MOMprof_StdDev=MOMprof_StdDev, iDmean=iDmean,
iSIPplot=iSIPplot, nr_SIPs_prof_mean=nr_SIPs_prof_mean, t_vec_GVplot=t_vec_GVplot_vec, ntGV_vec=ntGV_vec
)
#produce summary plot with final moment value
if (iFinal == 1):
PS.PlotMomentsFinal(MOMFinal)
#produce summary plot with crossing times
if (iTimesCross == 1):
PS.PlotMomentsTimeCross(TimesCross)
#<<<<<<<< end routine plotMultipleSimulationsMOM(filepath_input)
def plotMultipleSimulationGV(filepath_input, igetSIPdataonly=0, nz_max=0, dz_vec=0):
nr_sims=len(filepath_input)
print('call of plotMultipleSimulationGV')
print('nr_sims:',nr_sims)
nEK_sip_plot=np.zeros([nr_sims,nr_inst_max,nr_GVplot_max,nr_sip_max])
mEK_sip_plot=np.zeros([nr_sims,nr_inst_max,nr_GVplot_max,nr_sip_max])
nr_SIPs_plot=np.zeros([nr_sims,nr_inst_max,nr_GVplot_max],dtype='int')
#GCCif (addSIPplot > 0)
if (igetSIPdataonly == 1):
#nz_max = 1000
nr_SIPs_av=np.zeros([nr_sims,nr_GVplot_max])
nr_SIPs_prof_av=np.zeros([nr_sims,nr_GVplot_max,nz_max])
t_vec_GVplot_vec=np.zeros([nr_sims,nr_GVplot_max])
ntGV_vec=np.zeros([nr_sims],dtype='int')
#GCCendif /* (addSIPplot > 0) */
nr_inst_vec=np.zeros(nr_sims,dtype='int')
Vtot_vec=np.zeros(nr_sims)
skal_m_vec=np.zeros(nr_sims)
for i_sim in range(0,nr_sims):
fp=filepath_input[i_sim]
print('process simulation:',fp)
# read SIP data
nEK_sip_plot_sim, mEK_sip_plot_sim,zDummy, nr_SIPs_plot_sim, nr_SIPs_prof_sim, t_vec_GVplot, V, skal_m = IO.readSingleSimulationSIPdata(fp,nr_sip_max=nr_sip_max)
nr_inst,nr_GVplot,nr_SIPs = nEK_sip_plot_sim.shape
Vtot_vec[i_sim] = V
skal_m_vec[i_sim] = skal_m
print('nr_inst,nr_GVplot,nr_SIPs',nr_inst,nr_GVplot,nr_SIPs)
print(nr_SIPs_plot_sim.shape)
nr_inst_vec[i_sim] = nr_inst
if (nr_inst > nr_inst_max):
print ('increase nr_inst_max: ', nr_inst, nr_inst_max )
nr_SIPs_plot[i_sim,0:nr_inst,:nr_GVplot] = nr_SIPs_plot_sim
nEK_sip_plot[i_sim,0:nr_inst,:nr_GVplot,:nr_SIPs]= nEK_sip_plot_sim
mEK_sip_plot[i_sim,0:nr_inst,:nr_GVplot,:nr_SIPs]= mEK_sip_plot_sim
#GCCif (addSIPplot > 0)
if (igetSIPdataonly == 1):
nr_SIPs_av[i_sim,:nr_GVplot]=np.mean(nr_SIPs_plot_sim,axis=0)
prof_dim=nr_SIPs_prof_sim.shape # nr_inst,nr_GVplot, nz
nz = prof_dim[2]
#GCCif (addSIPplot == 1)
nr_SIPs_prof_av[i_sim,:nr_GVplot,:nz]=np.mean(nr_SIPs_prof_sim,axis=0)*100/dz_vec[i_sim]
#GCCendif /* (addSIPplot == 1) */
#GCCif (addSIPplot == 2)
nr_SIPs_prof_av[i_sim,:nr_GVplot,:nz]=np.mean(nr_SIPs_prof_sim,axis=0)
#GCCendif /* (addSIPplot == 2) */
ntGV_vec[i_sim] = nr_GVplot
t_vec_GVplot_vec[i_sim,:nr_GVplot]= t_vec_GVplot
#GCCendif /* (addSIPplot > 0) */
if (igetSIPdataonly == 0):
# Generate size distribution plot
for Times_elem in Times_list:
#GCCif (MOM_skal_m == 0)
PS.PlotGV(nEK_sip_plot,mEK_sip_plot,nr_SIPs_plot,t_vec_GVplot,
iMultipleSims=1,nr_inst_vec=nr_inst_vec,
label=label,title=title,ilastGVonly=ilastGVonly,iTimes_select=Times_elem,
V=Vtot_vec)
#GCCelse
PS.PlotGV(nEK_sip_plot,mEK_sip_plot,nr_SIPs_plot,t_vec_GVplot,
iMultipleSims=1,nr_inst_vec=nr_inst_vec,
label=label,title=title,ilastGVonly=ilastGVonly,iTimes_select=Times_elem,
V=Vtot_vec, skal_m=skal_m_vec)
#GCCendif /* else (MOM_skal_m == 0) */
else:
return nr_SIPs_av, nr_SIPs_prof_av, t_vec_GVplot_vec, ntGV_vec
#<<<<<<<< end routine plotMultipleSimulationGV(filepath_input)
#GCCif (DISCRETE >= 1)
def plotMultipleSimulationRStD(filepath_input,label=None,title=None):
nr_sims=len(filepath_input)
print('call of plotMultipleSimulationRStD')
print('nr_sims:',nr_sims)
mEK_sip_plot=np.zeros([nr_sims,nr_inst_max,nr_GVplot_max,nr_sip_max])
nr_inst_vec=np.zeros(nr_sims,dtype='int')
for i_sim in range(0,nr_sims):
fp=filepath_input[i_sim]
print('process simulation:',fp)
# read SIP data
nEK_sip_plot_sim, mEK_sip_plot_sim,zDummy, nr_SIPs_plot_sim, nr_SIPs_prof_sim, t_vec_GVplot, V, skal_m = IO.readSingleSimulationSIPdata(fp,nr_sip_max=nr_sip_max)
nr_inst,nr_GVplot = nr_SIPs_plot_sim.shape
nr_inst_vec[i_sim] = nr_inst
#nr_SIPs_plot[i_sim,0:nr_inst,:] = nr_SIPs_plot_sim
mEK_sip_plot[i_sim,0:nr_inst,:nr_GVplot,:]= mEK_sip_plot_sim
# Generate size distribution plot
PS.PlotRelStdDev(mEK_sip_plot,t_vec_GVplot,iMultipleSims=1,nr_inst_vec=nr_inst_vec,label=label,title=title)
#<<<<<<<< end routine plotMultipleSimulationRStD(filepath_input,label=None,title=None)
#GCCendif /* (DISCRETE >= 1) */
def plotMultipleSimulationFLUX(filepath_input):
nr_sims=len(filepath_input)
print('call of plotMultipleSimulationFLUX')
print('nr_sims:',nr_sims)
nt_max=80
#GCCif (STD == 1)
iplotStdDev = 1
FluxIn_STD= np.zeros([1,nr_sims,nt_max,5])
FluxOut_STD = np.zeros([1,nr_sims,nt_max,5])
FluxInAcc_STD = np.zeros([1,nr_sims,nt_max,5])
FluxOutAcc_STD = np.zeros([1,nr_sims,nt_max,5])
#GCCendif /* (STD == 1) */
#GCCif (STD == 2)
iplotStdDev = 2
FluxIn_STD= np.zeros([2,nr_sims,nt_max,5])
FluxOut_STD = np.zeros([2,nr_sims,nt_max,5])
FluxInAcc_STD = np.zeros([2,nr_sims,nt_max,5])
FluxOutAcc_STD = np.zeros([2,nr_sims,nt_max,5])
#GCCendif /* (STD == 2) */
#GCCif (STD == 0)
iplotStdDev = 0
FluxIn_STD= None
FluxOut_STD = None
FluxInAcc_STD = None
FluxOutAcc_STD = None
#GCCendif /* (STD == 0) */
FluxIn= np.zeros([nr_sims,nt_max,5])
FluxOut = np.zeros([nr_sims,nt_max,5])
FluxInAcc = np.zeros([nr_sims,nt_max,5])
FluxOutAcc = np.zeros([nr_sims,nt_max,5])
print('FluxIn.shape,FluxOut.shape,FluxInAcc.shape,FluxOutAcc.shape')
print(FluxIn.shape,FluxOut.shape,FluxInAcc.shape,FluxOutAcc.shape)
nt_vec =np.zeros(nr_sims,dtype='int')
t_vec_Flux_vec=np.zeros([nr_sims,nt_max])
#read flux data
for i in range(nr_sims):
fp_in=filepath_input[i]
print('process simulation:',fp_in)
FluxIn_tmp,FluxOut_tmp,FluxInAcc_tmp,FluxOutAcc_tmp,t_vec_GVplot=IO.readSingleSimulationFluxes(fp_in)
nt_vec[i] = t_vec_GVplot.size-1
t_vec_Flux_vec[i,:nt_vec[i]] = t_vec_GVplot[1:]
FluxIn[i,:nt_vec[i],:]=np.mean(FluxIn_tmp,axis=0)
FluxOut[i,:nt_vec[i],:]=np.mean(FluxOut_tmp,axis=0)
FluxInAcc[i,:nt_vec[i],:]=np.mean(FluxInAcc_tmp,axis=0)
FluxOutAcc[i,:nt_vec[i],:]=np.mean(FluxOutAcc_tmp,axis=0)
#GCCif (STD == 1)
FluxIn_STD[0,i,:nt_vec[i],:]=np.std(FluxIn_tmp,axis=0)
FluxOut_STD[0,i,:nt_vec[i],:]=np.std(FluxOut_tmp,axis=0)
FluxInAcc_STD[0,i,:nt_vec[i],:]=np.std(FluxInAcc_tmp,axis=0)
FluxOutAcc_STD[0,i,:nt_vec[i],:]=np.std(FluxOutAcc_tmp,axis=0)
#GCCendif /* (STD == 1) */
#GCCif (STD == 2)
FluxIn_STD[:,i,:nt_vec[i],:]=np.abs(np.percentile(FluxIn_tmp,[10,90],axis=0)-np.squeeze(FluxIn[i,:nt_vec[i],:]))
FluxOut_STD[:,i,:nt_vec[i],:]=np.abs(np.percentile(FluxOut_tmp,[10,90],axis=0)-np.squeeze(FluxOut[i,:nt_vec[i],:]))
FluxInAcc_STD[:,i,:nt_vec[i],:]=np.abs(np.percentile(FluxInAcc_tmp,[10,90],axis=0)-np.squeeze(FluxInAcc[i,:nt_vec[i],:]))
FluxOutAcc_STD[:,i,:nt_vec[i],:]=np.abs(np.percentile(FluxOutAcc_tmp,[10,90],axis=0)-np.squeeze(FluxOutAcc[i,:nt_vec[i],:]))
#GCCendif /* (STD == 2) */
#GCCif (STD > 0)
print('FluxIn_STD.shape,FluxOut_STD.shape,FluxInAcc_STD.shape,FluxOutAcc_STD.shape')
print(FluxIn_STD.shape,FluxOut_STD.shape,FluxInAcc_STD.shape,FluxOutAcc_STD.shape)
#GCCendif /* (STD > 0) */
#generate flux plot
PS.PlotFluxesTime(FluxIn,FluxOut,FluxInAcc,FluxOutAcc,t_vec_Flux_vec,
iMean=1,iMultipleSims=1,io_sep=1,nt_vec=nt_vec,
label=label,title=title,text=text,
iplotStdDev = iplotStdDev, FluxIn_STD = FluxIn_STD, FluxOut_STD = FluxOut_STD, FluxInAcc_STD = FluxInAcc_STD, FluxOutAcc_STD = FluxOutAcc_STD)
def plotMultipleSimulationGVout(filepath_input,infalling=0):
nr_sims=len(filepath_input)
print('call of plotMultipleSimulationGVout')
print('nr_sims:',nr_sims)
nEK_sip_out=np.zeros([nr_sims,nr_inst_max,nr_sip_max])
mEK_sip_out=np.zeros([nr_sims,nr_inst_max,nr_sip_max])
tEK_sip_out=np.zeros([nr_sims,nr_inst_max,nr_sip_max])
nr_SIPs_out=np.zeros([nr_sims,nr_inst_max],dtype='int')
nr_inst_vec=np.zeros(nr_sims,dtype='int')
#A_vec=np.zeros(nr_sims)
for i_sim in range(0,nr_sims):
fp=filepath_input[i_sim]
print('process simulation:',fp)
# read SIP data
nEK_sip_out_sim, mEK_sip_out_sim, tEK_sip_out_sim, nr_SIPs_out_sim, Tsim, A = IO.readSingleSimulationSIP_outfalling_data(fp,infalling)
nr_inst,nr_SIPs = nEK_sip_out_sim.shape
#A_vec[i_sim] = A_sim
print('nr_SIPs_out_sim', nr_SIPs_out_sim)
print('nr_inst,nr_SIPs', nr_inst, nr_SIPs)
print(nr_SIPs_out_sim.shape)
nr_inst_vec[i_sim] = nr_inst
if (nr_inst > nr_inst_max):
print ('increase nr_inst_max: ', nr_inst, nr_inst_max )
nr_SIPs_out[i_sim,0:nr_inst] = nr_SIPs_out_sim
nEK_sip_out[i_sim,0:nr_inst,:nr_SIPs]= nEK_sip_out_sim/A
mEK_sip_out[i_sim,0:nr_inst,:nr_SIPs]= mEK_sip_out_sim
# Generate size distribution plot
PS.PlotGV(nEK_sip_out,mEK_sip_out,nr_SIPs_out,[3600],iMultipleSims=1,nr_inst_vec=nr_inst_vec,label=label,title=title, outfallingGV=1+infalling)
#<<<<<<<< end routine plotMultipleSimulationGVout(filepath_input)
#GCCendif /* (IMODE == 2) */
#--------------- end defintion of subroutines --------------------------------------------------------------------
#GCCif (IMODE == 1)
print('filepath',filepath)
#GCCif (GV_meanTD == 1)
plotSingleSimulationGV(filepath)
#GCCendif /*(GV_meanTD == 1)*/
#GCCif (rSIGMA_time >= 1)
plotSingleSimulationRStD(filepath)
#GCCendif /* (rSIGMA_time >= 1) */
#GCCif (MOM_meanTD >= 1)
plotSingleSimulationMOM(filepath)
#GCCendif /* (MOM_meanTD >= 1) */
#GCCif (COLUMN == 1)
#GCCif (FLUX_time == 1)
plotSingleSimulationFLUX(filepath)
#GCCendif /* (FLUX_time == 1) */
#GCCif (MOM_prof == 1)
plotSingleSimulationMOM(filepath,iprof=1,itot=0)
#GCCendif /* (MOM_prof == 1) */
#GCCif (GV_out == 1)
plotSingleSimulationGVout(filepath)
#GCCendif /*(GV_out == 1)*/
#GCCif (GV_in == 1)
plotSingleSimulationGVout(filepath, infalling=1)
#GCCendif /*(GV_in == 1)*/
#GCCendif /* (COLUMN == 1) */
#GCCendif /* (IMODE == 1) */
#GCCif (IMODE == 2)
#GCCif (GV_meanTD == 1)
plotMultipleSimulationGV(filepath_input)
#GCCendif /*(GV_meanTD == 1)*/
#GCCif (rSIGMA_time >= 1)
plotMultipleSimulationRStD(filepath_input)
#GCCendif /* (rSIGMA_time >= 1) */
#GCCif (MOM_meanTD >= 1 || MOM_prof == 1)
iProf=0
iTime=0
iFinal=0
iTimesCross=0
#GCCif (MOM_meanTD >= 1)
iTime=1
#GCCif (FINAL == 1)
iFinal=1
iTime=0
#GCCendif /* (FINAL == 1) */
#GCCif (TIMECROSS == 1)
iTimesCross=1
iTime=0
#GCCendif /* (TIMECROSS == 1) */
#GCCendif /* (MOM_meanTD >= 1) */
#GCCif (MOM_prof == 1)
iProf=1
#GCCendif /* (MOM_prof == 1) */
plotMultipleSimulationsMOM(filepath_input,iTime=iTime,iProf=iProf,iFinal=iFinal,iTimesCross=iTimesCross)
#GCCendif /* (MOM_meanTD >= 1 || MOM_prof == 1) */
#GCCif (FLUX_time == 1)
plotMultipleSimulationFLUX(filepath_input)
#GCCendif /* (FLUX_time == 1) */
#GCCif (GV_out == 1)
plotMultipleSimulationGVout(filepath_input)
#GCCendif /*(GV_out == 1)*/
#GCCif (GV_in == 1)
plotMultipleSimulationGVout(filepath_input, infalling=1)
#GCCendif /*(GV_in == 1)*/
#GCCendif /* (IMODE == 2 ) */
#print(os.path.dirname(os.path.abspath(__file__)))