PRO errorcalcNCAR
close,/all
tek_color

;******************************************************************************************
;IDL program to do error calculations using the real K matrix and detail
;files output from SFIT2
;
;Written by Rodica and Bec, with various pieces from or based on SFIT2ERS code by
;Steven Wood (distributed with SFIT2)
;
;Inputs that are currently manually entered include the name of the K file, detail file,
;top and bottom of partial column of interest and information for the determination of the1
;Sx file for calculating Sa (Sa determined by SFIT2 is at too higher precision (e-100!) for
;IDL, giving strange eigenvalues.

;Still to do:
; Loop so that multiple partial columns can be calculated simultaneously
; Useful output files
; Better Sx determination

;RB and RL, Jan 2007

;Updates:

; Added temperature calculation part (requires use of separate program doptpert to create
; the Ktemp matrix) (RB 25/1/08)
; Fixed partial column averaging kernel calculation (RB 29/1/08)
; Program now reads apriori uncertainty and correlation areas from detail file, and the
; our_Sa matrix now handles non-total-column correlation (RB 29/1/08)

; Major overhaul to work with NCAR files and Jim's perturbation programs, as well as
;setting up to read input file rather than having to hard-code. RB 07/2010
; Removed most of the plots from this code into a separate plotting function at the end
; And set it up to call the various functions if I haven't run them in advance

;*******************************************************************************************

;Read the required input parameters from input

junk=''
openr,1,'input'
readf, 1, junk ;header
readf,1,junk ;directory
blah=strsplit(junk,/extract)
dir=blah(0)
readf,1,junk ;gas
blah=strsplit(junk,/extract)
gas=blah(0)
readf,1,junk; site
blah=strsplit(junk,/extract)
site=blah(0)
readf,1,junk; top of column
blah=strsplit(junk,/extract)
toplayer=double(blah(0))
readf,1,junk; bottom of column
blah=strsplit(junk,/extract)
bottomlayer=double(blah(0))
readf,1,junk; snr
blah=strsplit(junk,/extract)
realSNR=double(blah(0))
readf,1,junk; line intensity
blah=strsplit(junk,/extract)
lineintensityuncert=double(blah(0))
readf,1,junk; line width
blah=strsplit(junk,/extract)
linewidthuncert=double(blah(0))
readf,1,junk; SZA
blah=strsplit(junk,/extract)
SZAuncert=double(blah(0))
readf,1,junk; ILS modulation
blah=strsplit(junk,/extract)
ILSuncert=double(blah(0))
readf,1,junk; ILS phase
blah=strsplit(junk,/extract)
phaseuncert=double(blah(0))
close, 1

; only edit these if you need to change the names

Kfil=dir+'/K.out.def' ; The name of the K file - put path in too, if not in same directory
Afil=dir+'/AK.out.def' ; The name of the AK file - put path in too, if not in same directory
Dfil=dir+'/cdetail.def' ; as per K but for detail file
Ktempfilename=dir+'/Kb_temp.cov'; the name of the Ktemp file, if you're going to do a temperature error
Klineintensityfilename=dir+'/Kb_linestrength.cov'; the name of the Klinestrength file, if you're going to do a line parameter estimate
Klinewidthfilename=dir+'/Kb_linewidth.cov'; the name of the Kwidth file, if you're going to do a line parameter estimate
KSZAfilename=dir+'/Kb_pointing.cov'; the name of the KSZA file, if you're going to do a SZA estimate
KILSfilename=dir+'/Kb_ILS.cov'; the name of the KILS file, if you're going to do a ILS estimate
Kphasefilename=dir+'/Kb_phase.cov'; the name of the Kphase file, if you're going to do a phase error estimate

if site eq 'TAB' then fastlayers=43; the number of layers in the fastcode file
if site eq 'MLO' then fastlayers=40; the number of layers in the fastcode file
if site eq 'FL0' then fastlayers=42; the number of layers in the fastcode file


if bottomlayer eq 0 then begin
 if site eq 'TAB' then bottomlayer=0.23 ;these may need to be changed as need rounded z from "z" read in
 if site eq 'MLO' then bottomlayer=3.40
 if site eq 'FL0' then bottomlayer=1.61
 if site eq 'EUR' then bottomlayer=0.61
endif

;************************
dum=''; a dummy string variable, used frequently but always overwritten!
Stemp=0.0
Slineintens=0.0
Slinewidth=0.0
SSZA=0.0
SILS=0.0
Sphase=0.0

;Start by reading inportant information from the K file

openr,1,Kfil
readf,1,dum
readf,1,m,n ; gives size of the matrices
n=fix(n)    ;force these to be integers SWW 11/04
m=long(m)  ; allow for very large spectra    SWW 11/04
K=dblarr(m,n)
readf,1,K
readf,1,dum
SeInv=dblarr(m)
readf,1,SeInv
readf,1,dum
SaInv=dblarr(n,n)
readf,1,SaInv
Sa=dblarr(n,n)
readf,1,dum
readf,1,Sa
;print,Sa
close,1
;

;Now do the calculations to determine A and D (D is also called G, or the gain matrix)
;
;Note that these A, D and K are relative to the apriori (so must be scaled), and that they
;include the components from the full state vector, not just the layers of interest (which will be
;identified later in the code).

KT=transpose(K)
KTSeInv=KT
for i=0,n-1 do begin
 for j=0L,m-1 do begin    ; force j to be long SWW 11/04
   KTSeInv(i,j)=KT(i,j)*SeInv(j)
 endfor
endfor
D=invert(KTSeInv#K+SaInv)#KTSeInv
A=D#K

;Now to get parameters from the detail file

openr,1,Dfil

; first get the number of layers - perhaps there is a better way of doing this, but I don't know it!

nlayers=0;
junk=''
dum=''
for i=0,12-1 do readf,1, junk
 readf,1,junk
 readf,1,nlayers
; print, nlayers
i=0
close, 1

; now open it again, and find the parameters we want!

openr,1,Dfil

;read the zpt information

aprcounter=0 ;a fix for the problem of 2 retrieved profiles
WHILE NOT EOF(1) DO BEGIN

readf, 1, junk

if (junk eq ' RELATIVE UNCERTAINTIES OF THE A PRIORI PROFILE' and aprcounter eq 0) then begin ;fix to stop it reading in the second set of a priori uncertainties
; print, 'Found apriori uncertainties'
 apr_uncert=fltarr(nlayers)
 readf,1, apr_uncert
 readf,1, dum

 if (dum eq '') then begin
  zflg=0 ;flag for no correlation
  print, 'no correlation'
 endif else begin
  temp1=strsplit(dum,/extract); take out the half width
  if temp1(6) eq 'GAUSSIAN' then zflg=1
  if temp1(6) eq 'EXPONENTIAL' then zflg=2
  half_width=fix(temp1(10), type=4)
  readf,1,dum
  temp1=strsplit(dum,/extract); now take out min altitude
  min_corr_alt=fix(temp1(7), type=4)
  readf,1,dum
  temp1=strsplit(dum,/extract); now take out max altitude
  max_corr_alt=fix(temp1(7), type=4)
  ;print, 'max_corr_alt=', max_corr_alt
 ;print, half_width
 endelse
aprcounter=aprcounter+1
endif


if (junk eq ' PRESSURE-TEMPERATURE PROFILE FROM TAPE9') then begin
; print, 'i found pressure-temperature profile from tape9'
 readf,1,dum
 zpt=fltarr(3, nlayers)
 readf,1, zpt
endif

;get the air mass data

if (junk eq' MASS PATH DATA FROM TAPE7 IN CM*ATM' ) then begin
 ;print,'i found massp data'

 mass=fltarr(nlayers)
 blaharr=fltarr(nlayers)
 atc=fltarr(nlayers); total column AK
 sum=fltarr(nlayers)

 for i=0,3 do readf,1, dum ;read through the junk
  readf,1,blaharr
  readf,1,junk
  readf,1,blaharr
  readf,1,junk
  readf,1, mass

endif


;get the signal to noise ratio

if (junk eq ' DEFAULT SIGNAL-TO-NOISE' ) then begin
 ;print,'i found snr'
 readf,1,dum
 ;print,dum
 readf,1,dum
 ;print,dum
 junkarr=fltarr(4,1)
 readf,1, junkarr
 snr=junkarr(1)
 ;print,'snr=', snr
endif


;identify which part of each matrix corresponds to what, by attaching an identifying parameter

if (junk eq ' PRINT OUT OF PARAMETERS FOR EACH ITERATION:') then begin

  readf,1, junk
  nlines=0
  nlines=n/5
  parameters=strarr(n)

  for i=0, nlines-1 do begin
    readf,1, junk
    temp=strsplit(junk, /extract)
    parameters(i*5+0)=temp(0)
    parameters(i*5+1)=temp(1)
    parameters(i*5+2)=temp(2)
    parameters(i*5+3)=temp(3)
    parameters(i*5+4)=temp(4)
  endfor

  numleft=n-(5*nlines)
  readf,1, junk
  temp=strsplit(junk, /extract)

  for p=0, numleft-1 do begin
   parameters(i*5+p)=temp(0+p)
  endfor

  print, 'parameters are: ',parameters

endif

;now get the apriori and retrieved profile (and all these other things, all saved in an array called "retrieval"
 retrieval=fltarr(7, nlayers)

if (junk eq '  Z[km] ZBAR[km] APRIORI_VMR RETRIEVED_VMR  SIGMA_VMR  APRIORI_COL  RETRIEVE_COL') then begin

  header=''
  junk=header
  readf,1, retrieval
  endif

endwhile



close, 1


;***********************************************************************************
;Now to actually use some of this stuff!  First, identify which are the correct
;bits of the A and D matrices for our gas
;note newA and newD are going to be used quite a bit. They are still not scaled at this point

il=where(parameters eq gas,nil)

newA=fltarr(nlayers, nlayers)
newK=fltarr(m,nlayers)
newD=fltarr(nlayers,m)
newSaInv=fltarr(nlayers,nlayers)

newA=A[il(0):il(nil-1), il(0):il(nil-1)]
;print,'newA= ',newA
newK=K[0:(m-1),il(0):il(nil-1)]
newD=D[il(0):il(nil-1),0:(m-1)]
newSainv=SaInv[il(0):il(nil-1),il(0):il(nil-1)]
Sav=Sa[il(0):il(nil-1), il(0):il(nil-1)]
;note: at Toronto we found that IDL can't cope with the level of precision in the K file, and Sav contains some
;very odd values. Seems to be OK with UNIX

;Other important things:

;apriori profile:
apr=retrieval(2,*)
;print,'apriori ', apr

;z profile
z=zpt(0,*)

;now scale A, so that it is no longer the internal relative to apriori

Avmr=fltarr(nlayers,nlayers)

for i=0, nlayers-1 do begin
 for j=0, nlayers-1 do begin
 Avmr(i,j)= newA(i,j)*(apr(i)/apr(j))
  endfor
endfor
;print,'Avmr=',Avmr

;And use this, with the airmass info, to determine the total column averaging kernel

for j=0, nlayers-1 do begin
 for i=0, nlayers-1 do begin
  sum(i)= mass(i)*Avmr(i,j)
 endfor
atc(j)=(total(sum))*1./mass(j)
endfor

;print,'atc=', atc


;***************************************
;Error calculations

;From Rogers:
;Sm=D#Sm#DT
;Ss=(A-I)#Sa#(A-I)T
;
;***************************************

;create an identity matrix
Id=identity(nlayers,/double)

;print,'Id=',Id

AvI=fltarr(nlayers,nlayers)
AvI=Avmr-Id
AvIT=transpose(AvI)


;**** Ss ****

;A should be scaled. All of these are "relative" in that A, I and our_Sa are all dimensionless

Ss=AvI#Sav#AvI  ;note this is equivalent to (AVI##our_Sa##AVIT) ;previously used our_Sa
;print,'Ss=',Ss


;**** now for the measurement error ******

Sm=fltarr(nlayers,nlayers)
DT=Transpose(newD)
;print,'DT=',DT

;create the Se matrix, from our signal to noise ratio (assuming non-correlated errors,
; therefore off-diagonal elements are zero)

Se=fltarr(m,m)
dif=fltarr(m)+1./realSNR ;set to use SNR from input file, rather than one read from cdetail

for i=0,m-1 do begin
  Se(i,i)=dif(i)^2
endfor
;print,Se= ',Se

;**** Sm ****

Sm=(newD)#Se#DT
;print,'Sm=',Sm

;Note Sm and Ss are complete matrices

;output these covariance matrices
formatcode='('+string(nlayers)+'E20.10)'
SSfilename=dir+'SS.covariance'
openw,1,SSfilename
printf,1,'Smoothing error SS=(A-I)#Sa#(A-I)T'
for i=0, nlayers-1 do begin
 printf,1,SS(*,i), format=formatcode
endfor
close,1

SMfilename=dir+'SM.covariance'
openw,1,SMfilename
printf,1,'Measurement error SM=G#Se#GT'
for i=0, nlayers-1 do begin
 printf,1,sm(*,i), format=formatcode
endfor
close,1

;;*********************************************************

;Now for some of the partial column calculations

;detrmine the part of the column that we're interested in

topindexpart=where(z le toplayer)
topindex=min(topindexpart)

bottomindexpart=where(z le bottomlayer)
bottomindex=min(bottomindexpart)

print, 'top=', z(topindex)
print, 'bottom= ', z(bottomindex)

;create an array with ones where we're interested, and zeros where we arent.
partialrange=fltarr(nlayers)

for x=0, nlayers-1 do begin
 if (x ge topindex and x le bottomindex) then begin
  partialrange(x)=1.0
 endif else begin
  partialrange(x)=0.0
 endelse
endfor

;and a version with zeros where we're not interested, and the a priori where we are

newapr=apr*partialrange

;determining the partial column amount from the retrieved profile

nlay=1 ;at the moment, only one top and bottom layer

ret_pcol=fltarr(nlay)
for i=0, nlay-1 do begin
 ret_pcol(i)=total(retrieval(3,topindex:bottomindex)*mass(topindex:bottomindex))
endfor


;now determine the partial column averaging kernel, using the scaled Avmr
;then determine the partial column kernel by incorporating the density

sum_pcol=fltarr(nlayers)
ak_pcol=fltarr(nlayers)

for j=0, nlayers-1 do begin
 for i=0, nlayers-1 do begin
  sum_pcol(i)= partialrange(i)*mass(i)*Avmr(i,j)
 endfor
 ak_pcol(j)=(total(sum_pcol))*1./mass(j)
endfor

;and calculate the DOFS for the partial column
pcol_dofs=trace(newA[topindex:bottomindex, topindex:bottomindex])


;Bec's plot of individual layer columns, profiles and partial columns

set_plot,'ps'
device,filename='AK_plot.ps'
device,/color
!p.multi=[0,2,2]
!p.font=1

plot, Avmr(0,*), z, xrange=[-0.5,1.5], yrange=[0,120], tit='Averaging kernels (scaled)'
for i=0, nlayers-1 do begin
 oplot, Avmr(i,*), z, color=i+2
 if (i mod 2 eq 0) then xyouts,0.6,z(i),string(format='(i3)',z(i)),color=i+2, charsize=0.5
 if (i mod 2 eq 1) then xyouts,0.8,z(i),string(format='(i3)',z(i)),color=i+2, charsize=0.5
endfor
sens=fltarr(nlayers)
for i=0, nlayers-1 do begin
 sens(i)=total(Avmr(*,i))
endfor
oplot, sens, z, linestyle=2


dofs=trace(newA)
dofs1=strtrim(string(dofs),1)
dofs2=strmid(dofs1,0,5)

plot, newA(0,*), z, xrange=[-0.5,1.5], tit='Averaging kernels (unscaled)', xtitle='dotted=sum'
 xyouts,-0.4,85,'DOFS='
 xyouts,-0.1,85,dofs2
for i=0, nlayers-1 do begin
 oplot, newA(i,*), z, color=i+2
 if (i mod 2 eq 0) then xyouts,0.6,z(i),string(format='(i3)',z(i)),color=i+2, charsize=0.5
 if (i mod 2 eq 1) then xyouts,0.8,z(i),string(format='(i3)',z(i)),color=i+2, charsize=0.5
endfor
sens=fltarr(nlayers)
for i=0, nlayers-1 do begin
 sens(i)=total(newA(*,i))
endfor
oplot, sens, z, linestyle=2


 plot, retrieval(3,*), z, title='Profiles'
 oplot, retrieval(2,*), z, color=2
 xyouts,2e-6,85,'a priori', color=2
 xyouts,2e-6,79,'retrieved'


plot, ak_pcol(*), z, tit='Column averaging kernels', xrange=[-0.2,1.8], yrange=[0,120]
oplot, atc, z, color=2

dofs1=strtrim(string(dofs),1)
dofs2=strmid(dofs1,0,5)

pcol_dofs1=strtrim(string(pcol_dofs),1)
pcol_dofs2=strmid(pcol_dofs1,0,5)

zb=strtrim(string(z(0)),1)
zb1=strmid(zb,0,5)

zt=strtrim(string(z(nlayers-1)),1)
zt1=strmid(zt,0,5)

zbi=strtrim(string(z(bottomindex)),1)
zbid=strmid(zbi,0,5)

zti=strtrim(string(z(topindex)),1)
ztid=strmid(zti,0,5)

 xyouts,0.5,80,string('From: ', zt1, '  to   ', zb1), charsize=0.8,color=2
 xyouts,-0.4,80, string('DOFS=', dofs2), charsize=0.8, color=2
 xyouts,0.5,70,string('From: ', zbid, '  to  ', ztid), charsize=0.8;,color=i+2
 xyouts,-0.4,70, string('DOFS=', pcol_dofs2), charsize=0.8;, color=i+2

device,/close
;----------------------------------------------------

;*** Now calculate the partial column errors ****
; this gets a bit trickier

;partial column = apriori * air mass

Ss_err=fltarr(nlay)
Sm_err=fltarr(nlay)
St_err=fltarr(nlay)

; *** For the smoothing error ***

; there are two ways to do this - either using the already scaled Avmr as the A, in which case
; the partial column averaging kernel needs to be multiplied with the airmass, or using the unscaled
; A (done here), in which case the partial column averaging kernel needs a factor of a priori*airmass.

; Strictly, the averaging kernel needs a normalization to make it dimensionless. The Sx matrix for
; the smoothing error calculation then needs to be multiplied into partial column units. It turns out
; that not normalizing the averaging kernel leaves the partial column units in the A-I part, and if you
; don't scale the Sx, the numbers come out the same. But for the purpose of making the calculations
; correct, the following method has the normalizations in them too.


;1st, create two vectors that'll be used for scaling

; partdens: the density profile in the area we're interested in, zeros elsewhere
partdens=newapr*mass
; densprf: the density profile over the whole column, utilizing the a priori
densprf=apr*mass

;2nd, do the scaling that makes everything real units instead of relative to the apriori

Avmr_pcol=transpose((partdens)#newA) ; scaling A by the density*airmass where we're interested
Ivmr_pcol=partialrange*densprf; scaling I (in this case 1s where we're interested, 0s elsewhere) the same way

;Atest=Avmr_pcol
;Itest=Ivmr_pcol

;and normalize (2nd half of conversion of A(&I) to density/column or vmr)

for j=0,nlayers-1 do begin
   Avmr_pcol(j)=Avmr_pcol(j)/densprf(j)
   Ivmr_pcol(j)=Ivmr_pcol(j)/densprf(j)
endfor

; Calculate A-I
A_I_pcol=Avmr_pcol-Ivmr_pcol

;Now to get the units right, Sa also needs to be multiplied by the density!

Sa_scaled=dblarr(nlayers, nlayers) ;does this make sense at all?
for i=0, nlayers-1 do begin
 for j=0, nlayers-1 do begin
 Sa_scaled(i,j)=Sav(i,j)*densprf(i)*densprf(j)
 endfor
endfor


;Finally, determine the smoothing error

Ss_err=sqrt((A_I_pcol)#Sa_scaled#(A_I_pcol))
;print, 'ss_err=', ss_err

;Ss_err2=sqrt((Atest-Itest)#Sav#(Atest-Itest)) ; this confirms the alternate method listed above
;print, 'ss_err2=', ss_err2
;stop

; *** For the measurement error ***

;As for A, we need to make D relative to the density, by multiplying by partdense

D_pcol=transpose(partdens#newD);transpose puts it the right way around for calculating sm

Sm_err=sqrt(D_pcol#Se#D_pcol)


; *** And the total Sm + SS error ***

St_err=sqrt((Sm_err)^2+(Ss_err)^2)

;;print it all out
;
print, 'Column amount=',ret_pcol
print, 'Sm(%)=', (Sm_err/ret_pcol)*100
print, 'Ss(%)=', (Ss_err/ret_pcol)*100
;print, 'St(%)=', (St_err/ret_pcol)*100
;
print, 'DOFS (partial column) =', pcol_dofs
print, 'DOFS (total column) = ', trace(newA)

;******************************************************************

;Piece of code to look at "interference errors"
;slightly strange because there are two separate parts of the interference error in this - those before the vmr AK
;- these are the retrieval parameters like wavenumber shift, and those after - these are the interfering species.
; calculate each separately:

;Interference error 1 - retrieval parameters
;remembering that n=number of paramters in A and il=where(parameters eq gas,nil)

;get the right bit of the matrix (remembering IDL reads down, then across)
int1_AK=fltarr(nlayers, il(0))
int1_AK=A[il(0):il(nil-1), 0:il(0)-1]

;make the Se matrix

Se_int1=fltarr(il(0), il(0))

for i=0, il(0)-1 do begin

if (parameters(i) eq 'BckGrdSlp') then Se_int1(i,i)=0.1
if (parameters(i) eq 'BckGrdCur') then Se_int1(i,i)=0.1
if (parameters(i) eq 'IWNumShft') then Se_int1(i,i)=0.1
if (parameters(i) eq 'SWNumShft') then Se_int1(i,i)=0.1
if (parameters(i) eq 'EmpApdFcn') then Se_int1(i,i)=0.2
if (parameters(i) eq 'SolLnShft') then Se_int1(i,i)=0.5
if (parameters(i) eq 'SPhsErr') then Se_int1(i,i)=0.2

endfor

;Calculate Interference errors part 1 to get the error covariance matrix for the profile

Sint1=int1_AK#Se_int1#transpose(int1_AK)

;diagonal elements of Sint1
diagSint1=fltarr(nlayers)
Sint1sd=fltarr(nlayers)
for i=0,nlayers-1 do begin
  diagSint1(i)=Sint1(i,i)
  Sint1sd(i)=Sqrt(Sint1(i,i))
 endfor


;then for the partial column, multiplying the averaging kernel by density in the form of partialrange to get it right for the partial column

int1_AK_pcol=partdens#int1_AK

;calculating the error due to interference errors part 1

S_int1_err=sqrt(int1_AK_pcol#Se_int1#transpose(int1_AK_pcol))

print,'Sint1 (retrieval params) (%) = ', (S_int1_err/ret_pcol)*100

;the covariance matrix of this

Sint1filename=dir+'Sint1.covariance'
openw,1,Sint1filename
printf,1,'Interfering errors, including background slope, wavenumbershift, solar lines etc Sint1=AK_int1#Se_int#(AK_int1)T'
for i=0, nlayers-1 do begin
 printf,1,Sint1(*,i), format=formatcode
endfor
close,1


;Interference error 2 - interfering species
int2_AK=fltarr(nlayers, ((n-1)-il(nil-1)))

int2_AK=A[il(0):il(nil-1), (il(nil-1)+1):n-1]

Se_int2=fltarr(((n-1)-il(nil-1)), ((n-1)-il(nil-1)))

for i=0, ((n-1)-il(nil-1))-1 do begin
 Se_int2(i,i)=1.
endfor

Sint2=int2_AK#Se_int2#transpose(int2_AK)

;diagonal elements of Sint2
diagSint2=fltarr(nlayers)
Sint2sd=fltarr(nlayers)
for i=0,nlayers-1 do begin
  diagSint2(i)=Sint2(i,i)
  Sint2sd(i)=Sqrt(Sint2(i,i))
 endfor

Sint1filename=dir+'Sint2.covariance'
openw,1,Sint1filename
printf,1,'Interfering errors due to interfering species (assuming 100% uncert in interfering species) Sint2=AK_int2#Se_int#(AK_int2)T'
for i=0, nlayers-1 do begin
 printf,1,Sint2(*,i)
endfor
close,1

;then for the partial column, multiplying the averaging kernel by density in the form of partialrange to get it right for the partial column

int2_AK_pcol=partdens#int2_AK

;calculating the error due to interference errors part 2

S_int2_err=sqrt(int2_AK_pcol#Se_int2#transpose(int2_AK_pcol))

print,'Sint2 (intf. spec.)(%) = ', (S_int2_err/ret_pcol)*100


;
;********************************************************
;A small appendum, for doing the error due to temperature calculation. Note that this requires that
;you have already run dtemplbl and created the K_temp matrix, which is saved in one giant list in
;Ktemp_cols.txt

num=0
print, 'I can now continue to an error due to temperature calculation.'
print, 'Select from the following options:'
print, '1. Do not do error due to temperature calculation'
print, '2. I have run dtemplbl and Kb_temp exists. Do calculation'
print, '3. I have not run dtemplbl. Run dtemplbl then do calculation'
read,num, prompt='Enter selection: '

case num of
1: GOTO, line_intensity_calculation
2: GOTO, temp_continue
3: dtemplbl, fastlayers=fastlayers, nlayers=nlayers
else: begin
 print, 'Not a valid option. Exiting'
 stop
 end
endcase

temp_continue:

;Build the Ktemp matrix from program to read in all the columns and put them in the matrix

Ktemp=fltarr(fastlayers, m)
dum=''

openr,1,Ktempfilename
readf,1,Ktemp
Ktemp=transpose(Ktemp)
close,1

;make the Sb matrix

;Estimate the uncertainty in the temperature profile at each altitude

;uncertainty: Above 50km, 9K based on top NCEP height error, All uncertainties below based on
; NCEP error estimates: 40-50km = 7K, 35-40km =6K, 30-35km=5K, 2K below due to radisondes **Note, NCEP temp estimates
;have ~17-30km at 2.5-3km; this might be more realistic***

case site of
'EUR': begin
   li=[9, 9, 9, 9, 9, 9, 7, 7, 7, 7, 6, 6, 5, 5, 2, 2, 2, 2, 2, 2,$
   2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
 end
'TAB': begin
   li=[9,9,9,9,9,9,9,9,7,7,7,6,6,5,5,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,$
   2,2,2,2,2,2,2,2,2,2,2,2,2]
 end
'MLO': begin
   li=[9,9,9,9,9,9,9,9,9,9,9,9,7,7,7,7,7,6,6,6,5,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]
 end
 'FL0': begin
   li=[9,9,9,9,9,9,9,9,7,7,7,6,6,5,5,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,$
   2,2,2,2,2,2,2,2,2,2,2,2,2]
 end
else: print,'Do not have the temperature Sb yet. Edit errorcalcncar.pro'
endcase

;Make the matrix
Sb=fltarr(fastlayers,fastlayers)

;Formulation suggested by Brian Connor, with a correlation between layers. The following is from Becs thesis work.
;rho=0.0
;zc=8 ;for 8km correlated

;for i=0,nlayers-1 do begin
; for j=0, nlayers-1 do begin
;  rho=exp(-(abs(z(j)-z(i))/zc))
;  Sb(i,j)=rho*li(i)*li(j)
; endfor
;endfor

;Simple version
for i=0,fastlayers-1 do begin
  Sb(i,i)=li(i)*li(i)
endfor

;For D, use D_pcol - this is for the partial column. For full covariance see below

;Now to calculate Stemp column covariance
Stemp=sqrt(D_pcol#Ktemp#Sb#transpose(Ktemp)#D_pcol) ;no transpose on final D_pcol cos of way IDL handles vectors in matrix mult

;and the Stemp full covariance
Stempcov=newD#Ktemp#Sb#transpose(Ktemp)#transpose(newD)

print, 'Stemp(%)=', (Stemp/ret_pcol)*100

;output the covariance matrix

STempfilename=dir+'Stemp.covariance'
openw,1,STempfilename
printf,1,'Forward model temperature error Stemp=G#Kb#Sb#Kb#G'
for i=0, nlayers-1 do begin
 printf,1,Stempcov(*,i), format=formatcode
endfor
close,1

;********************************************************
line_intensity_calculation:

num=0
print, 'I can now continue to error calculations due to uncertainty in the line intensity and air broadened half width.'
print, 'Select from the following options:'
print, '1. Do not do error calculation due to spectral parameters'
print, '2. I have run dspec and Kb_linestrength and Kb_line width already exist. Do calculation'
print, '3. I have not run dspec. Run dspec then do calculation'
read,num, prompt='Enter selection: '

case num of
1: GOTO, SZA_calculation
2: GOTO, line_continue
3: dspec
else: begin
 print, 'Not a valid option. Exiting'
 stop
 end
endcase

line_continue:

;Build the K matrix from program to read in all the columns and put them in the matrix

Klineintensity=fltarr(m)
dum=''
openr,1,Klineintensityfilename
readf,1,Klineintensity
close,1

;make the Sb matrix - easy with only one parameter - it's a single number!

SbLI=(lineintensityuncert)^2

;Now calculate Slineintens
;column
Slineintens=sqrt(D_pcol##transpose(Klineintensity)##SbLI##Klineintensity##transpose(D_pcol))
;covariance
Slineintensitycov=newD#Klineintensity#SbLI#transpose(Klineintensity)#transpose(newD)

print, 'Uncertainty due to line intensity (Slineintens) = ', Slineintens
print, 'Slineintens(%)=', (Slineintens/ret_pcol)*100

;output the covariance matrices

Slineintensfilename=dir+'Slineintens.covariance'
openw,1,Slineintensfilename
printf,1,'Line intensity error S lineintens=G#Kb#Sb#Kb#G'
for i=0, nlayers-1 do begin
 printf,1,Slineintensitycov(*,i), format=formatcode
endfor
close,1

;line_strength_calculation:

;Build the K matrix from program to read in all the columns and put them in the matrix

Klinewidth=fltarr(m)
dum=''
openr,1,Klinewidthfilename
readf,1,Klinewidth
close,1

;make the Sb matrix - easy with only one parameter - it's a single number!

SbLW=(linewidthuncert)^2

;Now calculate Slinewidth
;column
Slinewidth=sqrt(D_pcol##transpose(Klinewidth)##SbLW##Klinewidth##transpose(D_pcol))

;and the Slinewidth covariance
Slinewidthcov=newD#Klinewidth#SbLW#transpose(Klinewidth)#transpose(newD)

print, 'Uncertainty due to line width (Slinewidth) = ', Slinewidth
print, 'Slinewidth(%)=', (Slinewidth/ret_pcol)*100

;output the covariance matrices

Slinewidthfilename=dir+'Slinewidth.covariance'
openw,1,Slinewidthfilename
printf,1,'Line width error S linewidth=G#Kb#Sb#Kb#G'
for i=0, nlayers-1 do begin
 printf,1,Slinewidthcov(*,i), format=formatcode
endfor
close,1


;***************************************************************
SZA_calculation:

num=0
print, 'I can now continue to error calculations due to uncertainty in the SZA/pointing.'
print, 'Select from the following options:'
print, '1. Do not do error calculation due to SZA'
print, '2. I have run dpoint and Kb_pointing already exists. Do calculation'
print, '3. I have not run dpoint. Run dpoint then do calculation'
read,num, prompt='Enter selection: '

case num of
1: GOTO, ILS_calculation
2: GOTO, sza_continue
3: dpoint
else: begin
 print, 'Not a valid option. Exiting'
 stop
 end
endcase

sza_continue:

;Build the K matrix from program to read in all the columns and put them in the matrix

KSZA=fltarr(m)
dum=''
openr,1,KSZAfilename
readf,1,KSZA
close,1

;make the Sb matrix - easy with only one parameter - it's a single number!

SbSZA=(SZAuncert)^2


;Now to actually calculate Ssza
;column

SSZA=sqrt(D_pcol##transpose(KSZA)##SbSZA##KSZA##transpose(D_pcol))

;and the Ssza covariance

Sszacov=newD#KSZA#SbSZA#transpose(KSZA)#transpose(newD)


;print, 'Uncertainty due to SZA (SSZA) = ', SSZA
print, 'SSZA(%)=', (SSZA/ret_pcol)*100

;output the covariance matrices

SSZAfilename=dir+'Spointing.covariance'
openw,1,SSZAfilename
printf,1,'Forward model SZA/pointing error S pointing=G#Kb#Sb#Kb#G'
for i=0, nlayers-1 do begin
 printf,1,Sszacov(*,i), format=formatcode
endfor
close,1

;****************************************************************

ILS_calculation:

num=0
print, 'I can now continue to error calculations due to uncertainty in the ILS modulation efficiency.'
print, 'Select from the following options:'
print, '1. Do not do error calculation due to ILS modulation efficiency'
print, '2. I have run dlshap and Kb_ILS already exists. Do calculation'
print, '3. I have not run dlshap. Run dlshap then do calculation (note reply will be needed)'
read,num, prompt='Enter selection: '

case num of
1: GOTO, phase_calculation
2: GOTO, ils_continue
3: dlshap
else: begin
 print, 'Not a valid option. Exiting'
 stop
 end
endcase

ils_continue:

;Build the K matrix from program to read in all the columns and put them in the matrix

KILS=fltarr(m)
dum=''
openr,1,KILSfilename
readf,1,KILS
close,1

;make the Sb matrix - easy with only one parameter - it's a single number!

SbILS=(ILSuncert)^2

;Now to actually calculate SILS
;column
SILS=sqrt(D_pcol##transpose(KILS)##SbILS##KILS##transpose(D_pcol)) ;there are some oddities to get the matrix multiplication to work. I'm concerned.

;and the Sils covariance

Silscov=newD#KILS#SbILS#transpose(KILS)#transpose(newD)


print, 'Uncertainty due to ILS mod efficiency (SILS) = ', SILS
print, 'SILS(%)=', (SILS/ret_pcol)*100

;output these covariance matrices

SILSfilename=dir+'SILS.covariance'
openw,1,SILSfilename
printf,1,'Forward model ILS modulation efficiency S_ils=G#Kb#Sb#Kb#G'
for i=0, nlayers-1 do begin
 printf,1,Silscov(*,i), format=formatcode
endfor
close,1

;****************************************************************

phase_calculation:

num=0
print, 'I can now continue to error calculations due to uncertainty in the ILS phase error.'
print, 'Select from the following options:'
print, '1. Do not do error calculation due to ILS phase error'
print, '2. I have run dphase and Kb_phase already exists. Do calculation'
print, '3. I have not run dphase. Run dphase then do calculation (note reply will be needed)'
read,num, prompt='Enter selection: '

case num of
1: GOTO, summary
2: GOTO, phase_continue
3: dphase
else: begin
 print, 'Not a valid option. Exiting'
 stop
 end
endcase

phase_continue:

;Build the K matrix from program to read in all the columns and put them in the matrix

Kphase=fltarr(m)
dum=''
openr,1,Kphasefilename
readf,1,Kphase
close,1

;make the Sb matrix - easy with only one parameter - it's a single number!

Sbphase=(phaseuncert)^2

;Now to calculate Sphase
;column
Sphase=sqrt(D_pcol##transpose(Kphase)##Sbphase##Kphase##transpose(D_pcol)) ;there are some oddities to get the matrix multiplication to work.

;and the Sphase covariance
Sphasecov=newD#Kphase#Sbphase#transpose(Kphase)#transpose(newD)


print, 'Uncertainty due to ILS phase error (Sphase) = ', Sphase
print, 'Sphase(%)=', (Sphase/ret_pcol)*100

;output the covariance matrix

Sphasefilename=dir+'Sphase.covariance'
openw,1,Sphasefilename
printf,1,'Forward model ILS phase error S_phase=G#Kb#Sb#Kb#G'
for i=0, nlayers-1 do begin
 printf,1,Sphasecov(*,i), format=formatcode
endfor
close,1

;****************************************************************

summary:

print, ''
print, 'SUMMARY:'
print, 'Sm(%)=', (Sm_err/ret_pcol)*100
print, 'Ss(%)=', (Ss_err/ret_pcol)*100
print, 'Stemp(%)=', (Stemp/ret_pcol)*100
print, 'Slineintens(%)=', (Slineintens/ret_pcol)*100
print, 'Slinewidth(%)=', (Slinewidth/ret_pcol)*100
print, 'Sint1 (retrieval params) (%) = ', (S_int1_err/ret_pcol)*100
print, 'Sint2 (intf. spec.) (%) = ', (S_int2_err/ret_pcol)*100
print, 'SSZA(%)=', (SSZA/ret_pcol)*100
print, 'SILS-mod(%)=', (SILS/ret_pcol)*100
print, 'SILS-phase(%)=', (Sphase/ret_pcol)*100

St_err1=sqrt((Sm_err)^2+(Stemp)^2+(S_int1_err)^2+(S_int2_err)^2+(SSZA)^2+(SILS)^2+(Sphase)^2)
print, 'Stot_random(%)=', (St_err1/ret_pcol)*100

St_err3=sqrt((Slineintens)^2+(Slinewidth)^2)
print, 'Stot_spectral(%)=', (St_err3/ret_pcol)*100

St_err2=sqrt((Sm_err)^2+(Stemp)^2+(S_int1_err)^2+(S_int2_err)^2+(SSZA)^2+(SILS)^2+(Sphase)^2+(Slineintens)^2+(Slinewidth)^2)
print, 'Stot_random+spectral(%)=', (St_err2/ret_pcol)*100


St_err=sqrt((Sm_err)^2+(Ss_err)^2+(Stemp)^2+(S_int1_err)^2+(S_int2_err)^2+(Slineintens)^2+(Slinewidth)^2+(SSZA)^2+(SILS)^2+(Sphase)^2)


print, 'Stot_incl_Ss(%)=', (St_err/ret_pcol)*100

print, 'DOFS (partial column) =', pcol_dofs
print, 'DOFS (total column) = ', trace(newA)

;put this into a text file too:

openw,1,'Error_summary'
printf,1,'Summary of column errors, '+string(bottomlayer)+' - '+string(toplayer)
printf,1,'Retrieved partial column = ', ret_pcol
printf,1, 'DOFS (partial column) =', pcol_dofs
printf,1, 'DOFS (total column) = ', trace(newA)
printf,1, ' '
printf,1, 'Sm(%)=', (Sm_err/ret_pcol)*100
printf,1, 'Ss(%)=', (Ss_err/ret_pcol)*100
printf,1, 'Stemp(%)=', (Stemp/ret_pcol)*100
printf,1, 'Slineintens(%)=', (Slineintens/ret_pcol)*100
printf,1, 'Slinewidth(%)=', (Slinewidth/ret_pcol)*100
printf,1, 'Sint1 (retrieval params) (%) = ', (S_int1_err/ret_pcol)*100
printf,1, 'Sint2 (intf. spec.) (%) = ', (S_int2_err/ret_pcol)*100
printf,1, 'SSZA(%)=', (SSZA/ret_pcol)*100
printf,1, 'SILS-mod(%)=', (SILS/ret_pcol)*100
printf,1, 'SILS-phase(%)=', (Sphase/ret_pcol)*100
printf,1, 'Stotal_random(%)=', (St_err1/ret_pcol)*100
printf,1, 'Stotal_spectral(%)=', (St_err3/ret_pcol)*100
printf,1, 'Stotal_random+spectral(%)=', (St_err2/ret_pcol)*100
printf,1, 'Stot_incl_Ss(%)=', (St_err/ret_pcol)*100
close,1


read, dum, prompt='Would you like to make plots (y or n)? '
if dum eq 'y' then begin
 errorplot, nlayers=nlayers, z=z, sa=sav, apr=retrieval(2,*)
endif

read, dum, prompt='Make error matrices (y or n)?'
if dum eq 'y' then begin

Srandom=fltarr(nlayers, nlayers)
Srandom=sqrt(Sm^2+Sint1^2+Sint2^2+Stempcov^2+Sszacov^2+SILScov^2+Sphasecov^2)
Srandomfilename=dir+'Srandom.covariance'
openw,1,Srandomfilename
printf,1,'Random error, measurement+interference+temperature+pointing+ILS covariance matrices, added in quadrature'
for i=0, nlayers-1 do begin
 printf,1,Srandom(*,i), format=formatcode
endfor
close,1

Ssystematic=fltarr(nlayers, nlayers)
Ssystematic=sqrt(Slineintensitycov^2+Slinewidthcov^2)
Ssystematicfilename=dir+'Ssystematic.covariance'
openw,1,Ssystematicfilename
printf,1,'Systematic error, spectral line intensity and air broadened half width'
for i=0, nlayers-1 do begin
 printf,1,Ssystematic(*,i), format=formatcode
endfor
close,1
endif

STOP
end