pro errorplot, nlayers=nlayers, z=z, sa=sa, apr=apr

IF( NOT KEYWORD_SET( z )) THEN BEGIN
    print, 'format = errorplot, nlayers=nlayers (int), z=z (fltarr[nlayers])'
    stop
ENDIF
IF( NOT KEYWORD_SET( nlayers )) THEN BEGIN
    READ, nlayers, PROMPT=' Enter number of retrieval layers (TAB=47): '
ENDIF

;***************************************
; Rather than have lots of error plotting as part of "errorcalc"  I've pulled it out and put it in here.
; Plots include error covariances etc
; RB July 2010
;***************************************
print, 'Making covariance error plots'

junk=''
;Read in the covariance matrices
Sm=fltarr(nlayers,nlayers)
Smout=fltarr(nlayers, nlayers)

Ss=fltarr(nlayers,nlayers)
Ssout=fltarr(nlayers, nlayers)

Sint1=fltarr(nlayers,nlayers)
Sint2=fltarr(nlayers,nlayers)


Stemp=fltarr(nlayers,nlayers)
Ssza=fltarr(nlayers,nlayers)
Slintens=fltarr(nlayers,nlayers)
Slwidth=fltarr(nlayers,nlayers)
Sphase=fltarr(nlayers,nlayers)
Sils=fltarr(nlayers,nlayers)

openr,1,'SM.covariance'
readf,1,junk
readf,1,sm
close,1

openr,1,'SM.out'
readf,1,junk
readf,1,smout
close,1

openr,1,'SS.covariance'
readf,1,junk
readf,1,ss
close,1

openr,1,'SS.out'
readf,1,junk
readf,1,ssout
close,1

openr,1,'Sint1.covariance'
readf,1,junk
readf,1,sint1
close,1

openr,1,'Sint2.covariance'
readf,1,junk
readf,1,sint2
close,1


openr,1,'SILS.covariance'
readf,1,junk
readf,1,sils
close,1

openr,1,'Stemp.covariance'
readf,1,junk
readf,1,stemp
close,1

openr,1,'Sphase.covariance'
readf,1,junk
readf,1,sphase
close,1

openr,1,'Slineintens.covariance'
readf,1,junk
readf,1,slintens
close,1

openr,1,'Slinewidth.covariance'
readf,1,junk
readf,1,slwidth
close,1

openr,1,'Spointing.covariance'
readf,1,junk
readf,1,ssza
close,1


Sm_diag=fltarr(nlayers)
for i=0,nlayers-1 do Sm_diag(i)=Sm(i,i)

Smout_diag=fltarr(nlayers)
for i=0,nlayers-1 do Smout_diag(i)=Smout(i,i)

Ss_diag=fltarr(nlayers)
for i=0,nlayers-1 do Ss_diag(i)=Ss(i,i)

Ssout_diag=fltarr(nlayers)
for i=0,nlayers-1 do Ssout_diag(i)=Ssout(i,i)

Sint1_diag=fltarr(nlayers)
for i=0,nlayers-1 do Sint1_diag(i)=Sint1(i,i)
Sint2_diag=fltarr(nlayers)
for i=0,nlayers-1 do Sint2_diag(i)=Sint2(i,i)
Stemp_diag=fltarr(nlayers)
for i=0,nlayers-1 do Stemp_diag(i)=Stemp(i,i)
Ssza_diag=fltarr(nlayers)
for i=0,nlayers-1 do Ssza_diag(i)=Ssza(i,i)
Sphase_diag=fltarr(nlayers)
for i=0,nlayers-1 do Sphase_diag(i)=Sphase(i,i)
Sils_diag=fltarr(nlayers)
for i=0,nlayers-1 do Sils_diag(i)=Sils(i,i)
Slintens_diag=fltarr(nlayers)
for i=0,nlayers-1 do Slintens_diag(i)=Slintens(i,i)
Slwidth_diag=fltarr(nlayers)
for i=0,nlayers-1 do Slwidth_diag(i)=Slwidth(i,i)

loadct,39, ncolors=12
set_plot,'ps'
device, filename='Error_plots.ps',/color
!p.multi=[0,2,1]

plot, sqrt(Sm_diag), z, tit='Sm, blue=SFIT SM.out, black=SM.covariance'
oplot, sqrt(Smout_diag), z, color=2

plot, sqrt(Ss_diag), z, tit='Ss, blue=SFIT SS.out, black=SS.covariance'
oplot, sqrt(Ssout_diag), z, color=2


!p.multi=[0,1,1]

plot, sqrt(Ss_diag), z, tit='Diagonal components of the covariance matrices', ytit='Altitude (km)', xtit='Error contribution';, xrange=[0,0.1]
oplot, sqrt(Sm_diag),z, color=2
oplot, sqrt(Sint1_diag),z, color=3
oplot, sqrt(Sint2_diag),z, color=4
oplot, sqrt(Stemp_diag),z, color=5
oplot, sqrt(Ssza_diag),z, color=6
oplot, sqrt(Sphase_diag),z, color=7
oplot, sqrt(Sils_diag),z, color=8
oplot, sqrt(Slintens_diag),z, color=9
oplot, sqrt(Slwidth_diag),z, color=10

xyouts, 0.15, 100, 'Ss'
xyouts, 0.15, 95, 'Sm', color=2
xyouts, 0.15, 90, 'Sinterfering retrieval params', color=3
xyouts, 0.15, 85, 'Sinterfering species', color=4
xyouts, 0.15, 80, 'Stemp', color=5
xyouts, 0.15, 75, 'Ssza', color=6
xyouts, 0.15, 70, 'Sphase', color=7
xyouts, 0.15, 65, 'Smodeff', color=8
xyouts, 0.15, 60, 'Slineintens', color=9
xyouts, 0.15, 55, 'Slinewidth', color=10

;random errors
Srand=sqrt(Sm^2+Sint1^2+Sint2^2+Stemp^2+Ssza^2+Sphase^2+Sils^2)
Ssyst=sqrt(Slintens^2+Slwidth^2)
Stot=sqrt(Sm^2+Sint1^2+Sint2^2+Stemp^2+Ssza^2+Sphase^2+Sils^2+Slintens^2+Slwidth^2+Ss^2)
Stot_ss=sqrt(Sm^2+Sint1^2+Sint2^2+Stemp^2+Ssza^2+Sphase^2+Sils^2+Slintens^2+Slwidth^2)

Srand_diag=fltarr(nlayers)
for i=0,nlayers-1 do Srand_diag(i)=Srand(i,i)
Ssyst_diag=fltarr(nlayers)
for i=0,nlayers-1 do Ssyst_diag(i)=Ssyst(i,i)
Stot_diag=fltarr(nlayers)
for i=0,nlayers-1 do Stot_diag(i)=Stot(i,i)
Stot_ss_diag=fltarr(nlayers)
for i=0,nlayers-1 do Stot_ss_diag(i)=Stot_ss(i,i)
Sa_diag=fltarr(nlayers)
for i=0,nlayers-1 do Sa_diag(i)=Sa(i,i)

plot, sqrt(Stot_diag)*100, z, tit='Error contributions using combined covariance matrices', ytit='Altitude (km)', xtit='Standard deviation (%)';, xrange=[0,0.1]
oplot, sqrt(Stot_ss_diag)*100, z, color=2
oplot, sqrt(Srand_diag)*100, z, color=4
oplot, sqrt(Ssyst_diag)*100, z, color=6
oplot, sqrt(Ss_diag)*100, z, color=8
oplot, sqrt(Sa_diag)*100,z, color=10
xyouts, 0.15, 100, 'Stotal'
xyouts, 0.15, 95,  'Stotal-Ss', color=2
xyouts, 0.15, 90, 'Srandom', color=4
xyouts, 0.15, 85, 'Ssystematic', color=6
xyouts, 0.15, 80, 'Ss', color=8
xyouts, 0.15, 75, 'A priori', color=10

plot, apr, z, tit='Apriori with errors from random error covariance matrix', xtit='VMR', ytit='Altitude'
for i=0, nlayers-1 do begin
 errbar=fltarr(2,2)
 errbar(0,0)=(apr(i)-sqrt(srand_diag(i))*apr(i))
 errbar(0,1)=(apr(i)+sqrt(srand_diag(i))*apr(i))
 errbar(1,0)=z(i)
 errbar(1,1)=z(i)
 print, errbar(0,0), errbar(1,0), errbar(0,1), errbar(1,1)
 oplot,errbar(0,*), errbar(1,*), color=10, thick=2
endfor
device,/close



set_plot,'ps'
device,/color, filename='Eigenvectors.ps'
;Now to look at the Eigenvectors and Eigenvalues, i.e. the error patterns
;
SST=transpose(Ss)
evalSm=HQR(ELMHES(Sm),/DOUBLE)
evecSm=EIGENVEC(Sm,evalSm,RESIDUAL=residual1)
evalSs=HQR(ELMHES(Ss),/DOUBLE)
evecSs=EIGENVEC(Ss,evalSs,RESIDUAL=residual2)
;
!p.multi=[0,4,2]
plot,evecSm(*,0),z,ytit='Height [km]', xtit='Error Pattern1', tit='Sm eigenvectors'
plot,evecSm(*,1),z,xtit='Error Pattern2'
plot,evecSm(*,2),z,xtit='Error Pattern3'
plot,evecSm(*,3),z,xtit='Error Pattern4'
plot,evecSm(*,4),z,ytit='Height [km]',xtit='Error Pattern5'
plot,evecSm(*,5),z,xtit='Error Pattern6'
plot,evecSm(*,6),z,xtit='Error Pattern7'
plot,evecSm(*,7),z,xtit='Error Pattern8'

plot,evecSs(*,0),z,ytit='Height [km]', xtit='Error Pattern1', tit='Ss eigenvectors'
plot,evecSs(*,1),z,xtit='Error Pattern2'
plot,evecSs(*,2),z,xtit='Error Pattern3'
plot,evecSs(*,3),z,xtit='Error Pattern4'
plot,evecSs(*,4),z,ytit='Height [km]', xtit='Error Pattern5'
plot,evecSs(*,5),z,xtit='Error Pattern6'
plot,evecSs(*,6),z,xtit='Error Pattern7'
plot,evecSs(*,7),z,xtit='Error Pattern8'
device,/close


end