#Langley data analysis

#for the reference spectrum see RRRRRRRRRRRRRR, otherwise do steps 1 and 2


#1111111111111111111111111111111111111111111111111111111111111111111111111111111111
#Step 1: Load level 1 files, extract SO routines, convert counts, and write file for windoas

fname='pan_prepareosscript.py'
File=open(fname);ss=File.readlines();File.close()
for iss in ss:exec(iss)

##Edit here-------------
#fixed parameters (edit only once)
pannum=6
icfname="c:\\X_Data\\pandora\\pandora6\\Pandora6_cf.txt"  #instrument calibration file
iofname="c:\\X_Data\\pandora\\pandora6\\Pandora6_of.txt"  #instrument operation file
pfad='C:\\X_Data\\pandora\\pandora6\\'  #level1 data file directory
#select staring and ending day
dd1=20081106
dd2=20090822
##dd1=20090807
##dd2=20090810
##----------------------

#read instrument operation file
sys_status=deepcopy(op.sys_status_ini)
sys_status[0][3]=iofname
res,spec_pars,hst_pars,tc_pars=io.get_instrparams(sys_status)
#read Instrument Calibration file
res,cal_data=io.get_calparams(icfname)
#Build last element of <cal_data>
cal_data=io.get_lastcalparam(cal_data,spec_pars,hst_pars)

#get wavelength from calibration
npix=spec_pars[5]  #number of pixels
pix=arange(npix)+1 #pixels fro 1 to ...
pixs=xfus.scale_x(pix,[npix])   #scaled pixels
lo=polyval(spec_pars[15],pixs)  #wavelengths as from calibration
i477=(lo>476)&(lo<478)     #wavelength rabge 476 to 478nm
ilonoi=lo>400     #wavelengths above 400nm
isl=lo<295  #for simple stray light correction
sysn="Pandora"

#loop over days
dds=range(dd1,dd2+1)
for idd in dds:
    dd=str(idd)
    if idd<20090000:
        stat="GSFC"
    else:
        stat="LaRC"
    weiter=True
    if int(dd[4:6])>12:
        weiter=False
    if weiter:
        #read header of level 1 file (to get column assignment)
        fname=pfad+sysn+str(pannum)+"_"+stat+"_"+dd+"_lev1.txt"
        res,instrnum,colind,locpars=io.read_lev1header(fname)
        if res[-9:]=="not exist":
            weiter=False
    if (weiter)and(res<>"OK"):
        ##if there is not header, then read the header of another file, that has a header;
        ##having no header is a bug in Version 0.4
        res,instrnum,colind,locpars=io.read_lev1header(pfad+"Pandora6_LaRC_20090807_lev1.txt")
    if weiter:
        latitude=locpars[2]
        longitude=locpars[3]
        #read open hole direct sun routines
        fname=pfad+sysn+str(pannum)+"_"+stat+"_"+dd+"_lev1.txt"
        tlc="SO"  
        res,a1=io.get_lev1data(fname,colind,tlc,ind="ALL")
        nmeas=len(a1)

        cca=[]  #this list will have the corrected count rates for each measurement
        noia=[]  #this list will have the uncertainties for each measurement
        info=[] #this list will have additional info for each measurement
        #info
        #  0   1   2   3   4   5   6    7    8
        #zeit dur sza saz f477 it ncy tboard noi

        #loop over measurements
        for i in range(nmeas):
            a1i=a1[i]
            res,aa,comment=io.combine_data(a1i)
            #if measurement is OK (i.e. not saturated or interrupted), continue
            if (res=="OK")and(comment==[]):
                #apply dark and non-linearity correction and copnvert to count rates
                cc,ecc,ecci=coco.convert_rc(aa,cal_data,steps=[1,2,3])
                #additional info
                ibc=aa[2][:,2]==1
                zeit=aa[1][ibc,0]  #start time of measurement at open hole
                dur=aa[1][ibc,3]   #duration of measurement at open hole
                zeitc=zeit+0.5*dur #central time of measurement
                sza,saz=astro.time2zaaz(zeitc,latitude,longitude)  #solar position
                it=aa[2][ibc,0]    #integraiton time
                ncy=aa[2][ibc,1]   #number of cycles
                tboard=aa[3][ibc,1]  #board temperature
                if tboard==999:
                    ##in version 0.3 the board temperature was by mistake called detector temperature
                    ##and therefore in another column
                    tboard=aa[3][ibc,0]
                noi=ecc/cc*100     #percentage noise
                mnoii=noi[:,ilonoi].mean(axis=1)  #mean over visible pixels
                f477=cc[0,i477].mean()  #intensity at 477nm
                #simple stray ligh correction
                sl=cc[0,isl].mean()
                ccc=cc[0,:]-sl
                cca.append(ccc)
                noia.append(noi[0,:])
                infoi=[zeit[0],dur[0],sza[0],saz[0],f477,it[0],ncy[0],tboard[0],mnoii[0]]
                info.append(infoi)
        #convert lists to arrays
        info=array(info)
        cca=array(cca)
        noia=array(noia)
        nvalid=len(cca)
        print dd,nmeas,nvalid
    
        #write file for windoas
        if nvalid>0:
            fname=pfad+sysn+str(pannum)+"_"+stat+"_"+dd
            exten=".spe"
            res,ind,File=io.open_file(fname+exten,"w")
            loczeit=info[:,0]-5*3600
            timestr=xfus.zeit2timestr(loczeit[0],meth=3)
            date=timestr[6:]+"/"+timestr[4:6]+"/"+timestr[:4]
            tag=floor(loczeit/3600./24.)   #day
            stu=(loczeit/3600./24.-tag)*24 #hour
            hh=cca
            iu=hh<1
            hh[iu]=1
            for i in range(nvalid):
                ss=("%.2f" % info[i,2])+" "+date+" "+("%.5f" % stu[i])+" "
                for ipix in range(npix):
                    ss=ss+("%.1f" % hh[i,ipix])+" "
                ss=ss+"\n"
                ss=ss.encode("latin_1")
                File.write(ss)
            File.close()
            #dump <info>
            exten=".ddf"
            res,ind,File=io.open_file(fname+exten,"w")
            dump(info,File)
            File.close()
Beep(800,100)

#figures for last analyzed day
clf()
plot(lo,cca.T)

#integration times
clf()
plot(info[:,5],'.')

#noise
clf()
plot(info[:,-1])

#temp
clf()
plot(info[:,7])





#222222222222222222222222222222222222222222222222222222222222
#Step 2: after running windoas for each day, load windoas retrievals

infoa=zeros((0,9))   #has the info-array for all days
wwa=zeros((0,13))   #has the ww-array for all days
for idd in dds:
    dd=str(idd)
    if idd<20090000:
        stat="GSFC"
    else:
        stat="LaRC"
    weiter=True
    fname=pfad+sysn+str(pannum)+"_"+stat+"_"+dd
    if int(dd[4:6])>12:
        weiter=False
    if weiter:
        exten=".asc"
        res,ind,File=io.open_file(fname+exten,"r")
        if res=="OK":
            ss=File.readlines()
            File.close()
        else:
            weiter=False
    if weiter:
        #remove head lines
        ss=ss[10:]
        nss=len(ss)
        #read windoas results
        ww=zeros((nss,13))
        #ww
        # 0               1    2        3
        # Fractional-day SZA Best-shift	SZA
        # 4   5   6     7       8    9   10   11      12
        #RMS h2o uh2o o3t223 uo3t223 o4 uo4 no2t293 uno2t293
        for i in range(nss):
            hss=ss[i].split()
            for j in range(13):
                ww[i,j]=float(hss[j])
        #load <info>
        exten=".ddf"
        res,ind,File=io.open_file(fname+exten,"r")
        info=load(File)
        File.close()
        #synchronize data
        nmeas=len(info)
        loczeit=info[:,0]-5*3600
        tag=floor(loczeit/3600./24.)   #day
        stu=(loczeit/3600./24.-tag)*24 #hour
        dstu=10/3600.
        stuww=(ww[:,0]-floor(ww[:,0]))*24
        indww=zeros(nmeas)  #gets a one for the valid spectra
        for i in range(nmeas):
            dt=abs(stuww-stu[i])
            if dt.min()<dstu:
                indww[i]=1
        iww=indww>0
        #print dd,len(info[iww,0]),len(ww)
        #add info and ww to infoa and wwa
        infoa=vstack((infoa,info[iww,:]))
        wwa=vstack((wwa,ww))
        
        
#watch data
col=0  #any of the selected species
titel=["NO2","O2O2","O3","H2O","rms"]
einh=["1e15molc/cm2","1e40molc2/cm4","1e15molc/cm2","1e20molc/cm2",""]
colw=[11,9,7,5,4]
qw=[1/1e15,1,1/1e15,1,1]
qf=[op.gas[1][2]*op.unitcf[1][1]/1e15,op.gas[2][2]*op.unitcf[1][2]/1e40,op.gas[0][2]*op.unitcf[1][1]/1e15,op.gas[5][2]*op.unitcf[1][3]/1e20,0.01]

clf()
x=infoa[:,0]
y=wwa[:,colw[col]]*qw[col]
H=plot(x,y,'b.')
H=getp(gca(),'xticklabels');setp(H,size=fosi)
H=getp(gca(),'yticklabels');setp(H,size=fosi)
xlabel('TIME',fontsize=fosigr)
if (col>=0)and(col<4):
    ylabel('RELATIVE '+titel[col]+' SLANT COLUMN ['+einh[col]+"]",fontsize=fosigr)
else:
    ylabel('RMS OF RESIDUALS',fontsize=fosigr)
title(sysn+' '+str(pannum)+', '+titel[col],fontsize=fosigr)
grid(True)

#get effective NO2 height assuming std conditions
igas=xfus.get_index("NO2",op.gas,0)
gasinfo=op.gas[igas]
qgas=atmos.change_unit(gasinfo[2],gasinfo[3])
#get effective height
heff=0
for ii in gasinfo[4]:
    proftype=ii[0]
    blfact=ii[1]
    h,teff,surc=atmos.get_hteff(proftype,blfact)
    heff=heff+ii[2]*h
#get direct sun AMF
ndata=len(infoa)
amf=zeros(ndata)
for i in range(ndata):
    amf[i]=atmos.amf_direct(infoa[i,2],heff)
    
rms=wwa[:,4]
qr=wwa[:,11]/op.unitcf[1][1]

clf()
plot(amf,qr,'.')

clf()
plot(amf,rms,'.')




#make minimum amount Langley calibration on filtered data
groupanz=200
perc=0.02   #percentile used
amflim=5  #limiting amf
rmslim=0.005  #rms-limit

ig=(rms<rmslim)&(amf<amflim)
clf()
x=amf[ig]
y=qr[ig]
plot(x,y,'.')
grid()

nm=len(x)  #number of filtered data
ngroup=nm/groupanz  #number of groups

#sort in amf
ind=x.argsort()
x=x[ind];y=y[ind]
#get percentile in each group
p=zeros((ngroup,2))
for i in range(ngroup):
    xin=x[i*groupanz:(i+1)*groupanz]
    yin=y[i*groupanz:(i+1)*groupanz]
    p[i,0]=xin.mean()
    p[i,1]=xfus.percentile(yin,perc)
#fit straight line in percentiles, assume no error in x-axis    
res,s=xfus.linfit(p[:,0],p[:,1])
xx=array([0,amflim])
yy=s[0]*xx+s[1]

#figure
clf()
H=plot(x,y,'b.',p[:,0],p[:,1],'r.',xx,yy,'g',p[:,0],p[:,1],'r.')
H[0].set_label('RELATIVE NO2 SLANT COLUMNS')
setp(H[0],markersize=10)
H[1].set_label(str(perc*100)+'%ile LEVEL FOR '+str(ngroup)+' BINS')
setp(H[1],markersize=15)
H[2].set_label('SC='+("%.2f" % s[0])+'DU*AMF'+("%.2f" % s[1])+'DU')
setp(H[2],linewidth=4)
setp(H[3],markersize=15)
H=getp(gca(),'xticklabels');setp(H,size=fosi)
H=getp(gca(),'yticklabels');setp(H,size=fosi)
xlabel('DIRECT SUN AMF',fontsize=fosigr)
ylabel('RELATIVE NO2 SLANT COLUMNS [DU]',fontsize=fosigr)
title('MODIFIED LANGLEY CALIBRATION PANDORA '+str(pannum),fontsize=fosigr)
font=FontProperties(size=fosigr)
H=legend(loc=2,prop=font)
grid(True)

fname=pfad+sysn+str(pannum)+"_MODLANCAL.png"  #choose other format (bmp, jpg) if you want
savefig(fname,
        dpi=None,facecolor='w',edgecolor='w',
        orientation='landscape',transparent=False)


#add reference slant column and divide by AMF
qref=s[1]
qa=qr-qref  #absolute slant columns
qv=qa/amf  #absolute vertical columns

ig=rms<rmslim
clf()
plot(amf[ig],qr[ig],'.',amf[ig],qa[ig],'.',amf[ig],qv[ig],'.')

#plot vertical columns for each day
loczeit=infoa[:,0]-5*3600
tag=floor(loczeit/3600./24.)   #day
d=hstack((1,tag[1:]-tag[:-1]))
ind=where(d>0)
tags=tag[ind[0]]
ntags=len(tags)  #number of days
stu=(loczeit/3600./24.-tag)*24 #hour

for itag in range(ntags):
    iin=(tag==tags[itag])&(rms<rmslim)
    x=stu[iin]
    y=qv[iin]
    if sum(iin)>0:
        clf()
        H=plot(x,y,'.')
        setp(H[0],markersize=masi)
        H=getp(gca(),'xticklabels');setp(H,size=fosi)
        H=getp(gca(),'yticklabels');setp(H,size=fosi)
        xlabel('LOCAL HOUR',fontsize=fosigr)
        ylabel('NO2 TOTAL COLUMNS [DU]',fontsize=fosigr)
        timestr=xfus.zeit2timestr(tags[itag]*24*3600,meth=3)
        title('PANDORA '+str(pannum)+' DATA '+timestr,fontsize=fosigr)
        axis([5,19.4,0,2])
        fname=pfad+'NO2_'+sysn+str(pannum)+"_"+timestr+".png"
        savefig(fname,
                dpi=None,facecolor='w',edgecolor='w',
                orientation='landscape',transparent=False)

        
        
#save data in ascii file
fname=pfad+"NO2_"+sysn+str(pannum)
exten=".txt"
res,ind,File=io.open_file(fname+exten,"w")
ss="UT-Time TotalNO2[DU] fitting-rms"
for i in range(ndata):
    timestr=xfus.zeit2timestr(infoa[i,0],meth=4)
    ss=ss+"\n"+timestr+" "+("%.2f" % qv[i])+" "+("%.4f" % rms[i])
ss=ss.encode("latin_1")
File.write(ss)
File.close()













#RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR
#pick the reference spectrum and save it for windoas

#do step 1 for 20090810, then continue here

#f477 versus noise
clf()
plot(info[:,4],log10(info[:,-1]),'.')

#pick the measurement with smallest noise
mnoi=info[:,-1].min()
ind=where(info[:,-1]==mnoi)
ind=ind[0][0]

clf()
plot(lo,cca[ind,:])

#save reference spectrum
exten=".ref"
fname=pfad+sysn+str(pannum)+"_"+stat+"_"+dd
ss=""
hh=cca[ind,:]
iu=hh<1
hh[iu]=1
for i in range(npix):
    ss=ss+("%.3f" % lo[i])+" "+("%.2f" % hh[i])+"\n"
ss=ss[:-1]
res,ind,File=io.open_file(fname+exten,"w")
if res=="OK":
    ss=ss.encode("latin_1");File.writelines(ss);File.close()

#write dispersion and slit function files for windoas    

#dispersion
fname=pfad+sysn+str(pannum)+"_dispersion.clb"
ss=""
nlo=len(lo)
for i in range(nlo):
    ss=ss+("%.3f" % lo[i])+"\n"
ss=ss[:-1]
res,ind,File=io.open_file(fname,"w")
if res=="OK":
    ss=ss.encode("latin_1")
    File.writelines(ss)
    File.close()

#slit function
fname=pfad+sysn+str(pannum)+"_resolution.slf"
resol=0.52
ss=""
for i in range(nlo):
    ss=ss+("%.3f" % lo[i])+" "+("%.2f" % resol)+"\n"
ss=ss[:-1]
res,ind,File=io.open_file(fname,"w")
if res=="OK":
    ss=ss.encode("latin_1")
    File.writelines(ss)
    File.close()