; analysis_2obs
; program to examine effects of analysing 2 idealised observation,
; given a forecast background field and assumed error covariance functions
; This solves the analysis problem by explicitly calculating the inverse of 
; the relevant 2x2 matrix
; ASI version revised 10 April 2002
; ---------------------------------------------------------------------------

function wf, x
; define background field at point x
; In this version wf=0.0 ALWAYS 

w=0.0

return,w
end

function ocor, x1, x2, Lo
; define correlation between two observations according to their separation
; use a SOAR function with length scale Lo (Lo=0 => uncorrelated errors)

d=abs(x1-x2)

if Lo gt 0.0 then begin
  c=(1+d/Lo)*exp(-d/Lo)
end else begin
  c=0.0
end    

return,c
end

function fcor, x1, x2, Lf
; define background correlation between two points x1, x2
; use a SOAR function with length scale Lf (Lf=0 => uncorrelated errors)

d=abs(x1-x2)

if Lf gt 0.0 then begin
  c=(1+d/Lf)*exp(-d/Lf)
end else begin
  c=0.0
  if d eq 0.0 then c=1.0
end    

return,c
end

; -----------------------------------------------------------------------------
; MAIN program
; -----------------------------------------------------------------------------

print,"ANALYSIS_2OBS - Toy Analysis System"
print,"Constructs 1-D analysis from 2 idealised obs"
print,"with specified values and error characteristics"
print," "
print,"You may change the analysis parameters at the prompt, or.."

; set x values at which equations are solved
x0=0.
xn=10.
nx=100 ; actually have nx+1 points numbered 0 - nx
dx=(xn-x0)/nx
x=dx*findgen(nx+1)
f=fltarr(nx+1)
a=fltarr(nx+1)
a1=fltarr(nx+1)
a2=fltarr(nx+1)

; initial observation values and locations
x1=4.0
o1=wf(x1)+1.0
x2=6.0
o2=wf(x2)-1.0

; define covariance functions by sigma and correlation lengths
; forecast (background)
sigf=1.
Lf=2.0
; observation
sigo=1.
Lo=0.0

; Start main loop
; Set observation positions & values
loop: print,"Press ENTER to retain previous values"
print," "
print,'Current value o1 of ob#1',o1
onew=''
read,'New o1 (between -2 and +2): ',onew
if keyword_set(onew) then o1=float(onew)
print,'Current position x1 of ob#1',x1
xnew=''
read,'New x1 (between 0 and 10): ',xnew
if keyword_set(xnew) then x1=float(xnew)
print,'Current value o2 of ob#2',o2
onew=''
read,'New o2 (between -2 and +2): ',onew
if keyword_set(onew) then o2=float(onew)
print,'Current position x2 of ob#2',x2
xnew=''
read,'New x2 (between 0 and 10): ',xnew
if keyword_set(xnew) then x2=float(xnew)
print,'o1, x1, o2, x2: ',o1, x1, o2, x2
; Set errors
print,'Observation error sigo (used for both obs)',sigo
 snew=''
read,'New sigo: ',snew
if keyword_set(snew) then sigo=float(snew)
print,'The observation error correlation depends on the distance between obs'
print,'Current observation error correlation length Lo: ',Lo
cnew=''
read,'New Lo (put Lo = 0.0 for uncorrelated errors): ',cnew
if keyword_set(cnew) then Lo=float(cnew)
print,'Background error sigf (constant)',sigf
snew=''
read,'New sigf: ',snew
if keyword_set(snew) then sigf=float(snew)
print,'Background error correlation length Lf',Lf
cnew=''
read,'New Lf (put Lf = 0.0 for uncorrelated errors): ',cnew
if keyword_set(cnew) then Lf=float(cnew)
print,'sigo, Lo, sigf, Lf: ',sigo, Lo, sigf, Lf

; forecast fields at ob locations
f1=wf(x1)
f2=wf(x2)
; O-F values
d1=o1-f1
d2=o2-f2
; weighted correlation (for 2 obs)
C=sigf*sigf*fcor(x1,x2,Lf)+sigo*sigo*ocor(x1,x2,Lo)
C=C/(sigf*sigf+sigo*sigo)
W=(sigf*sigf+sigo*sigo)*(1.-C*C)
W=sigf*sigf/W
; weight for a single ob
W1=sigf*sigf/(sigf*sigf+sigo*sigo)

for i=0,nx do begin
  xi=x(i)
  f(i)=wf(xi)
  cf1=fcor(xi,x1,Lf)
  cf2=fcor(xi,x2,Lf)
  z=cf1*(d1-C*d2)+cf2*(d2-C*d1)
  a(i)=f(i)+W*z
  a1(i)=f(i)+cf1*W1*d1
  a2(i)=f(i)+cf2*W1*d2
end

; plot results

stitle=string(sigo,Lo,sigf,Lf,$
format='("sigo=",f4.2," Lo=",f4.2," sigf=",f4.2," Lf=",f4.2)')

plot,x,f,yrange=[-2.,2.],$
title='Toy analysis (dotted 1 ob, dashed 2 obs)',$
subtitle=stitle
oplot,[x1,x2],[o1,o2],psym=2
oplot,x,a1,linestyle=1
oplot,x,a2,linestyle=1
oplot,x,a,linestyle=2

psfile=''
print,'If you wish to keep this plot, give the name of a Postscript file'
read,'Postscript file: ',psfile
if keyword_set(psfile) then begin
  psfile=strtrim(psfile,2)
  ; find any '.', and derive basename of file
  i = strpos(psfile,'.')
  if i ne -1 then psfile=strmid(psfile,0,i)
  print, 'postscript file: ',psfile+'.ps'        
  plopen,fn=psfile
  !p.font=0
  
  plot,x,f,yrange=[-2.,2.],$
  title='Toy analysis (dotted 1 ob, dashed 2 obs)',$
  subtitle=stitle
  oplot,[x1,x2],[o1,o2],psym=2
  oplot,x,a1,linestyle=1
  oplot,x,a2,linestyle=1
  oplot,x,a,linestyle=2

  plclose
  !p.font=-1

end

again=''
read,'Do you want to try again? (y or n): ', again
if keyword_set(again) then begin
  again=strmid(again,0,1)
  if again EQ 'y' or again EQ 'Y' then begin
    print,' '
    goto, loop
  end  
end 
print,'Goodbye' 

end