io	=0	;0=X-window, 1=ps, 2=eps, 3=color ps
form	=1	;0=landscape, 1=portrait
char	=0.6	;character size 
fignr	='9a'  	;figure number
plotname='oxnard_f'	;filename
ref	=' (Aschwanden et al. 2002)'  ;label at bottom of Fig.	
unit	=0	;window number
angle	=150
niter	=1
;goto,stereo

;Density distribution model
nx	=512
ny	=nx
wid	=nx/10.
x	=findgen(nx)
y	=findgen(ny)
model	=fltarr(nx,ny)
ns	=5
x0	=nx*[0.1,0.3,0.5,0.6,0.65]
y0	=ny*[0.7,0.65,0.6,0.4,0.2]
a0	=1.*[0.5,0.8,1.0,0.85,0.75]
wx	=nx*[0.1,0.12,0.1,0.08,0.05]
wy	=nx*[0.1,0.12,0.1,0.09,0.08]
for is=0,ns-1 do begin 
 for j=0,ny-1 do model(*,j)=model(*,j)+a0(is)*exp(-((x-x0(is))^2/(2.*wx(is)^2)+(y(j)-y0(is))^2/(2.*wy(is)^2)))
endfor

;STEREO-A,B view 
STEREO:
stereoa =fltarr(nx)
for j=0,ny-1 do stereoa(*)=stereoa(*)+model(*,j)

stereob =fltarr(nx)
if (angle mod 180) eq 0 then stereob=stereoa
if (angle mod 180) ne 0 then begin
 modelb	=rot(model,angle)
 x_	=fltarr(nx)+x(ny/2)
 y_	=y
 rotate1,x_,y_,xb,yb,nx/2,ny/2,angle
 for i=0,nx-1 do stereob(i)=total(modelb(i,*))
endif

;DISPLAY
fig_open,io,form,char,fignr,plotname,unit 
!p.position=[0.1,0.7,0.45,0.95]
!p.title='Model density distribution n(x,z)'
!p.charsize=char
!x.range=[0,nx]
!y.range=[0,ny]
!x.style=1
!y.style=1
nlevel	=10
level	=max(model)*(findgen(nlevel)+1.)/float(nlevel)
contour,model,x,y,level=level
oplot,[1,1]*(nx/2),[0,ny]
xyouts,nx/2.+10,10,'A',size=char
oplot,xb,yb
xyouts,xb(0),yb(0),'B',size=char
!noeras=1

;3D-RECONSTRUCTION
for iter=0,niter-1 do begin
 recon	=fltarr(nx,ny)
 a	=(angle/180.)*!pi
 if (iter eq 0) then stereoz=stereob
 if (iter ge 1) then begin
  diff  =max(abs(recb-stereob),i0)
  diffpeak=recb(i0)-stereob(i0)
  sign	=diffpeak/abs(diffpeak)
  corr	=-(0.5*sign*diff*exp(-(x-x(i0))^2/(2.*wid^2))>0)
  stereoz=stereob+corr
 endif
 for i=0,nx-1 do begin
  ishift=float(i-nx/2.)/tan(a)
  recon(i,*)=stereoa(i)*shift(stereoz/total(stereoz),-ishift)
 endfor
 recon=recon*total(model)/total(recon)	;renormalization

 reca =fltarr(nx)
 recb =fltarr(nx)
 for j=0,ny-1 do reca(*)=reca(*)+recon(*,j)
 if (angle mod 180) eq 0 then recb=reca
 if (angle mod 180) ne 0 then begin
  reconb	=rot(recon,angle)
  for i=0,nx-1 do recb(i)=total(reconb(i,*))
 endif
 diffabs=total(abs(recon-model))/total(model)
 xyouts,0.6,0.9-0.03*iter,'Diff='+string(diffabs),/normal
;oplot,x,recb
;oplot,x,stereoz,linestyle=2
endfor

!p.position=[0.1,0.4,0.45,0.65]
!p.title='STEREO A LOS_integration  n(x)'
!y.range=[0,max(stereoa)]
plot,x,stereoa,thick=3
oplot,x,reca
oplot,[1,1]*(nx/2),!y.range
xyouts,nx/2.+10,10,'A',size=char

!p.position=[0.55,0.4,0.9,0.65]
!p.title='STEREO B LOS_integration  n(x_rot)'
!y.range=[0,max(stereob)]
plot,x,stereob,thick=3
oplot,[1,1]*(nx/2),!y.range
xyouts,nx/2.+10,10,'B',size=char

!p.position=[0.1,0.1,0.45,0.35]
!p.title='Reconstructed density distribution n(x,z)'
!p.charsize=char
!x.range=[0,nx]
!y.range=[0,ny]
!x.style=1
!y.style=1
contour,recon,x,y,level=level

!p.position=[0.55,0.1,0.9,0.35]
!p.title='Difference [recon-model]'
!p.charsize=char
!x.range=[0,nx]
!y.range=[0,ny]
!x.style=1
!y.style=1
contour,recon-model,x,y,level=level
contour,model-recon,x,y,level=level,c_linestyle=1

fig_close,io,fignr,ref
end