pro stereo_loop_incl,index1,data1,index2,data2,savefile,nsm
;+
; Project     : SECCHI/STEREO
;
; Name        : STEREO_LOOP_INCL()
;
; Category    : 3D Reconstruction, Stereoscopy, Data Analsis
;
; Explanation : From a traced loop in STEREO image 1, the inclination angle 
;		of the loop plane is determined by correlation with stripes 
;		in STEREO image 2,
;		according to the method of "dynamic steroscopy" described in 
;               "3D-Stereoscopic Analysis of Solar Active Region Loops:
;               I.SoHO/EIT Observations at Temperatures of 1.0-1.5 MK"
;               (Aschwanden et al. 1999, ApJ 517, 977-989)
;
; Syntax      : IDL> stereo_loop_incl,index1,data1,index2,data2,savefile
;
; Inputs      : index1  - structure containing descriptors of data1
;               data1   - image at time t1 (reference time)
;               index2  - structure containing descriptors of data2
;               data2   - image at time t2 (e.g.from following or previous day)
;               savefile- string filename (e.g. 'test_001.sav')
;		nsm	- number of pixels for smoothing/highpass-filtering
;
; Outputs     : output parameters are written in savefile
;
; History     : 16-Aug-2006, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;---------------------------------------------------------------;
;               CHECK INPUT                                     ;
;---------------------------------------------------------------;
print,'INPUT : index1,data1,index2,data2'
help,index1,data1,index2,data2
dim1    =size(data1)
dim2    =size(data2)
if (dim1(0) eq 0) or (dim2(0) eq 0) then stop,'No image input!'
zmax    =max(data1)
zmin    =min(data1)

;---------------------------------------------------------------;
;               INPUT PARAMETERS                                ;
;---------------------------------------------------------------;
restore,savefile        ;-->x_range,y_range
fov_x   =(x_range(1,0)-x_range(0,0))
fov_y   =(y_range(1,0)-y_range(0,0))
cdelt1_ =fltarr(3)  &cdelt1_(1)=index1.cdelt1  &cdelt1_(2)=index2.cdelt1
cdelt2_ =fltarr(3)  &cdelt2_(1)=index1.cdelt2  &cdelt2_(2)=index2.cdelt2
crval1_ =fltarr(3)  &crval1_(1)=index1.crval1  &crval1_(2)=index2.crval1
crval2_ =fltarr(3)  &crval2_(1)=index1.crval2  &crval2_(2)=index2.crval2
crpix1_ =fltarr(3)  &crpix1_(1)=index1.crpix1  &crpix1_(2)=index2.crpix1
crpix2_ =fltarr(3)  &crpix2_(1)=index1.crpix2  &crpix2_(2)=index2.crpix2
solar_r =fltarr(3)  &solar_r(1)=index1.solar_r &solar_r(2)=index2.solar_r
pos_    =fltarr(3)     &pos_(1)=index1.solar_p    &pos_(2)=index2.solar_p
crota2_ =fltarr(3)  &crota2_(1)=index1.crota2  &crota2_(2)=index2.crota2
blong_  =fltarr(3)   &blong_(1)=index1.solar_l0 &blong_(2)=index2.solar_l0
blat_   =fltarr(3)    &blat_(1)=index1.solar_b0  &blat_(2)=index2.solar_b0
dateobs_=strarr(3)  &dateobs_(1)=index1.date_obs &dateobs_(2)=index2.date_obs
pa      =pos_(1)
p0a     =crota2_(1)
l0a     =blong_(1)
b0a     =blat_(1)
rsun    =696.;Mm
hfoot   =2.5 ;Mm
rmin    =1.+(hfoot/rsun)

;---------------------------------------------------------------;
;		IMAGE FILTERING 				;
;---------------------------------------------------------------;
smooth_flatedge,float(data1),lowpass1,nsm
smooth_flatedge,float(data2),lowpass2,nsm
filter1=float(data1)-lowpass1
filter2=float(data2)-lowpass2
nw	=nsm*3

;---------------------------------------------------------------;
;		DISPLAY STEREO IMAGES				;
;---------------------------------------------------------------;
io	=0
ysize	=800
xsize	=(ysize*(fov_x/fov_y))
window,0,xsize=xsize,ysize=ysize,retain=2
ncol	=2
nrow	=2
nimages	=ncol*nrow
qgap	=0.1
grid	=5.
index_ 	=concat_struct(concat_struct(concat_struct(index1,index2),index1),index2)
data_	=[[[data1]],[[data2]],[[filter1]],[[filter2]]]
x_label ='EW [arcseconds]'
y_label ='NS [arcseconds]'
z_range =fltarr(2,nimages)
for i=0,1 do z_range(*,i)=[zmin,zmax]
for i=2,3 do z_range(*,i)=[-2,+4]*nsm
fig_multi_pos,x_range,y_range,ncol,nrow,qgap,p_position
fig_multi_tv,index_,data_,x_range,y_range,z_range,p_position,ncol,nrow,io,$
        x_label,y_label,grid
xyouts,0.1,0.97,dateobs_(1),/normal,size=1.2 
xyouts,0.6,0.97,dateobs_(2),/normal,size=1.2 

;---------------------------------------------------------------;
;                calculate azimuth and baseline 		;
;---------------------------------------------------------------;
ns	=n_elements(xs)
az0	=0.
th0	=0.
daz	=0.
dbase	=0.
x1	=xs(0)/solar_r(1)
y1	=ys(0)/solar_r(1)
x2	=xs(ns-1)/solar_r(1)
y2	=ys(ns-1)/solar_r(1)
helio_trans,     pa,p0a,l0a,b0a,x1,y1,rmin,l1,b1,rsun
helio_trans,     pa,p0a,l0a,b0a,x2,y2,rmin,l2,b2,rsun
if (abs(l1-l2) ge 180) then begin
 l20=l2
 if (l1 gt l20) then l2=l2+360.
 if (l1 le l20) then l2=l2-360.
endif
arctan,(l2-l1),(b2-b1),az_rad,az
;cos_b	=(!pi/180.)*(b1+b2)/2.			;incorrect
 cos_b	=cos( (!pi/180.)*(b1+b2)/2. )		;corrected by Costis Gontikakis, April 24, 2002
base   =sqrt(((l2-l1)*cos_b)^2+(b2-b1)^2)*2.*!pi*rsun/360. ;Mm
print,'azimuth     az   =',az  ,' deg'
print,'baseline    base =',base,' Mm'
  
;---------------------------------------------------------------;
;		stereoscopy with varying inclination angle	;
;---------------------------------------------------------------;
!p.position=p_position(*,3)
!x.range=x_range(*,3)
!y.range=y_range(*,3)
plot,[0,0],[0,0]

  ind5  =1+long((ns-3)*findgen(5)/4.)
  dpix	=cdelt1_(1)*0.7/rsun ;Mm
  thmax	=89.
  th1	=-thmax
  th2	=+thmax
  dth	=  1. ;deg
  nth	=long((th2-th1)/dth+1)
  coalen_=fltarr(nth)
  p      =pos_(2)
  p0     =crota2_(2)
  l0     =blong_(2)
  b0     =blat_(2)
  r0	 =solar_r(2)
  cdelt1 =cdelt1_(2)
  cdelt2 =cdelt2_(2)
  crpix1 =crpix1_(2)
  crpix2 =crpix2_(2)
  crval1 =crval1_(2)
  crval2 =crval2_(2)
  coalen_max=0.0
  th_opt =0
  th_	=th1+dth*findgen(nth)
  for ith=0,nth-1 do begin
   coalen=0.
   arrow=' '
   th	=th_(ith)
   x    =xs(ind5)/replicate(solar_r(1),ns)
   y    =ys(ind5)/replicate(solar_r(1),ns)
help,x,y
   helio_loopcoord2,pa,p0a,l0a,b0a,rsun,az,th,x(0:1),y(0:1),rmin,xl,yl,zl
   helio_trans,     pa,p0a,l0a,b0a,x(0),y(0),rmin,l1,b1,rsun
   helio_loopcoord, xl,yl,zl,rsun,az,th,l1,b1,rmin,l,b,r
   if (min(r) lt rmin) or (max(r) gt 1.5) then goto,next_th

   helio_loopcoord2,pa,p0a,l0a,b0a,rsun,az,th,x,y,rmin,xl,yl,zl
help,xl,yl,zl
   helio_loopcoord, xl,yl,zl,rsun,az,th,l1,b1,rmin,l,b,r
help,l,b,r
   helio_trans2,    p,p0,l0,b0,l,b,r,xx,yy,zz,rsun
help,xx,yy,zz
   loop_interpol,dpix,xx,yy,xi,yi
   ns	=n_elements(xi)
   ix	=crpix1-1.+(xi*r0-crval1)/cdelt1
   iy	=crpix2-1.+(yi*r0-crval2)/cdelt2
   array_curve,filter2,ix,iy,nw,ixw_,iyw_,zw_
   zw_s=smooth(zw_,3)

   dwj	=long(nsm/2)
   jw1	=3*dwj
   jw2	=nw-1-3*dwj 
   for j=jw1,jw2 do begin
    if (j eq jw1) then begin
     zw3	=fltarr(ns,3) 
     for jj=j-3*dwj,j-1*dwj do zw3(*,0)=zw_s(*,jj)
     for jj=j-1*dwj,j+1*dwj do zw3(*,1)=zw_s(*,jj)
     for jj=j+1*dwj,j+3*dwj do zw3(*,2)=zw_s(*,jj)
    endif
    if (j gt jw1) then begin
      zw3(*,0)=zw3(*,0)-zw_s(*,j-3*dwj)+zw_s(*,j-1*dwj)
      zw3(*,1)=zw3(*,1)-zw_s(*,j-1*dwj)+zw_s(*,j+1*dwj)
      zw3(*,2)=zw3(*,2)-zw_s(*,j+1*dwj)+zw_s(*,j+3*dwj)
    endif
    ipeak=where(zw3(*,1) gt (zw3(*,0)+zw3(*,2))/2.,npeak)
    if (npeak ge 1) then begin
     istart=where(ipeak-shift(ipeak,+1) ne +1)
     iend  =where(ipeak-shift(ipeak,-1) ne -1)
     ilen  =iend-istart+1
     qlen  =max(ilen,im)/float(ns) 
     if (qlen gt coalen) then begin 
      coalen=qlen 
      jm    =j 
      ip    =ipeak(istart(im):iend(im))
     endif
    endif
    if (coalen gt coalen_max) then begin
     coalen_max=coalen
     th_opt=th
     xs_opt=xi*r0
     ys_opt=yi*r0
     arrow='<--'
     if (qlen gt 0.5) then oplot,(ixw_(ip,jm)-(crpix1-1.))*cdelt1+crval1,$
           (iyw_(ip,jm)-(crpix2-1.))*cdelt2+crval2,color=200,thick=2
     xstripe1=(ixw_(*,0)   -(crpix1-1.))*cdelt1+crval1
     xstripe2=(ixw_(*,nw-1)-(crpix1-1.))*cdelt1+crval1
     ystripe1=(iyw_(*,0)   -(crpix2-1.))*cdelt2+crval2
     ystripe2=(iyw_(*,nw-1)-(crpix2-1.))*cdelt2+crval2
    endif
    if ((th mod 15) eq 0) and (abs(th) ne 90) then begin
     oplot,xi*r0,yi*r0,color=255
     xyouts,xi(ns/2)*r0,yi(ns/2)*r0,string(th,'(i3)')+'!Uo!N',size=0.5
    endif
   
   endfor
   coalen_(ith)=coalen
   next_th:
  if (coalen eq 0) then print,'th=',th,format='(a,i4)'
  if (coalen gt 0) then print,'th=',th,' r(1)=',r(1),max(r),'  ns=',ns,$
    '  q_coal >',coalen,arrow,format='(a,i4,a,2f6.3,a,i4,a,f6.3,a)'
  endfor

  oplot,xs_opt,ys_opt,thick=3
  oplot,[xstripe1,reverse(xstripe2),xstripe1(0)],[ystripe1,reverse(ystripe2),ystripe1(0)],linestyle=2

;---------------------------------------------------------------;
;	stereoscopic inclination angle comparison		;
;---------------------------------------------------------------;
theta_range:
clearplot
window,2,retain=2
!x.title='Inclination angle' 
!y.title='Maximum length of coaligned segment'
!y.range=0
plot,th_,coalen_(*)
;read,'Enter smoothing     nsm=',nsm_ang
nsm_ang=5
coal2	=smooth_edge(coalen_(*),nsm_ang,0)
oplot,th_,coal2,thick=2
oplot,th1*[1,1],!y.crange,linestyle=1
oplot,th2*[1,1],!y.crange,linestyle=1
q2	=max(coal2,ith2)  ;maximum
for i2=ith2+1,nth-1 do if coal2(i2) gt coal2(i2-1) then goto,valley2
valley2:
for i1=ith2-1,0,-1 do if coal2(i1) gt coal2(i1+1) then goto,valley1
valley1:
th	=th_(ith2)
oplot,th*[1,1],[0,q2]
print,'inclination angle th=',th,' + ',dth

;---------------------------------------------------------------;
;	save results 						;
;---------------------------------------------------------------;
save,filename=savefile,x_range,y_range,xs,ys,nw,l1,b1,base,dbase,az,daz,th,dth
print,'Parameters written in file = '+savefile

end