pro stereo_simul,date_obs,days,res,baseline,l1,b1,az,th,index1,data1,index2,data2
;+
; Project     : STEREOSCOPY 
;
; Name        : STEREO_SIMUL()
;
; Purpose     : Test stereoscopy code 
;
; Category    : 3D data analysis, stereoscopy
;
; Explanation : simulates a "stereo pair" of images DATA and structures INDEX
;		that mimic a static coronal loop with prescribed geometry and 
;		is projected according to the solar rotation rate
;
; Syntax      : IDL> stereo_simul,date_obs,days,res,baseline,l1,b1,az,th,
;			index1,data1,index2,data2
;
; Examples    : IDL> ssw_packages,/mjastereo
;		 date_obs='2000-01-01T00:00:00.000Z
;		 days    =[-4,2.]  ;time difference in days
;		 res     =2.5      ;pixel size [in arcseconds]
;		 baseline=100.     ;loop footpoint separation [in arcseconds]
;		 l1      =30.      ;heliograph. long. of one footpoint [deg]
;		 b1      =40.      ;heliograph. lat.  of one footpoint [deg]
;		 az      =50.      ;azimuth angle of loop baseline to EW direction
;		 th      =60.      ;inclination angle of loop plane to vertical 
;		 stereo_simul,date_obs,days,res,baseline,l1,b1,az,th,
;			     index1,data1,index2,data2
;		 stereo_loop,index1,data1,index2,data2
;
; Inputs      : date_obs	=Date of observation 
;		days		=day difference to reference time [date_obs] 
;		res		=resolution of image [arcseconds]
;		baseline	=footpoint separation in heliographic degrees
;		l1		=heliographic longitude of 1st loop footpoint
;		b1		=heliographic latitude  of 1st loop footpoint
;		az		=aziumuth angle of loop baseline to EW direction
;		th		=inclination angle of loop plane
;
; Opt. Inputs : None
;
; Outputs     : index1		=structure of 1st stereo image 
;  		data1		=image     of 1st stereo image 
;  		index2		=structure of 2nd stereo image 
;  		data2		=image     of 2nd stereo image 
;
; Opt. Outputs: None
;
; Keywords    : None
;
; Common      : None
;
; Restrictions: None
;
; Side effects: None
;
; History     : Version 1,  14-OCT-1998,  Markus J. Aschwanden.  Written
;		update documentation, 7-NOV-1998, MJA
;
; Contact     : aschwand@lmsal.com 
;-

;_____________________________TIME GRID______________________________
nd	=n_elements(days)
tgrid	=timegrid(anytim(date_obs,/ints),days=days,/ccsds)
index	=sswfits_struct(nd,/addfits)

;_____________________________STRUCTURES : INDEX_____________________
fov	=2000.		;arcsecs (full disk image)
naxis1	=long(fov/res)
naxis2	=long(fov/res)
cdelt1	=res
cdelt2	=res
crpix1	=(naxis1+1.)/2. 
crpix2	=(naxis2+1.)/2.
crval1	=0.		
crval2	=0.		
xcen	=crval1+cdelt1*((naxis1+1)/2.-crpix1)	;=0
ycen	=crval2+cdelt2*((naxis2+1)/2.-crpix2)	;=0
eph     =get_sun(tgrid)
r       =reform(eph(1,*))         ;solar radius in arcseconds
l0      =reform(eph(10,*))        ;heliographic longitude of solar disk
b0      =reform(eph(11,*))        ;heliographic latitude of solar disk
p       =reform(eph(12,*))        ;position angle of solar axis
crota2	=-p			  ;image is rotated to solar north

index.instrume='Simulation'
index.date_obs=tgrid
index.naxis1  =naxis1
index.naxis2  =naxis2
index.cdelt1  =cdelt1
index.cdelt2  =cdelt2
index.crval1  =crval1
index.crval2  =crval2
index.crpix1  =crpix1
index.crpix2  =crpix2
index.crota2  =crota2
index.xcen    =xcen
index.ycen    =ycen
index.solar_r =r
index.solar_l0=l0
index.solar_b0=b0
index.solar_p =p

;_____________________________SOLAR DISK IMAGE_______________________
data	=bytarr(naxis1,naxis2,nd)
image0	=lonarr(naxis1,naxis2)
xx	=crval1+cdelt1*(findgen(naxis1)-(crpix1-1))
yy	=crval2+cdelt2*(findgen(naxis2)-(crpix2-1))

r0      =index(0).solar_r
for j=0,naxis2-1 do begin
 if (abs(yy(j)) lt r0) then begin
  xr	=sqrt(r0^2-yy(j)^2)
  ind	=where((xx ge -xr) and (xx le xr))
  image0(ind,j)=long(50+50*cos(sqrt(xx(ind)^2+yy(j)^2)/r0))
 endif
endfor

window,0,xsize=naxis1,ysize=naxis2,retain=2
loadct,3
!p.position=[0,0,1,1]
!x.range=[xx(0),xx(naxis1-1)]
!y.range=[yy(0),yy(naxis2-1)]
plot,[0,0],[0,0]
tv,bytscl(image0,min=0,max=255)

;_____________________________ROTATING LOOP___________________________ 
for iday=0,nd-1 do begin
image	=image0
eph     =get_sun(tgrid(iday))
r0      =eph(1)         ;solar radius in arcseconds
l0      =eph(10)        ;heliographic longitude of solar disk
b0      =eph(11)        ;heliographic latitude of solar disk
p       =eph(12)        ;position angle of solar axis
p0	=-p		;rotation angle of image 
ir	=r0/cdelt1	;solar radius in pixels

phi	=!pi*findgen(181)/180
rsun	=696.
rloop	=(baseline/2.)*0.7	;in Mm
xl	=rloop*cos(phi)
yl	=phi*0.
zl	=rloop*sin(phi)
r1	=696.+2.5
helio_loopcoord,xl,yl,zl,rsun,az,th,l1,b1,r1,l,b,r
helio_trans2,p,p0,l0,b0,l,b,r,x,y,z,rsun
xr	=(x/696.)*r0
yr	=(y/696.)*r0
print,'xr=',minmax(xr)
print,'yr=',minmax(yr)

;oplot,xr,yr
crpix1	=index(iday).crpix1
crpix2	=index(iday).crpix2
crval1	=index(iday).crval1
crval2	=index(iday).crval2
cdelt1	=index(iday).cdelt1
cdelt2	=index(iday).cdelt2
ii	=long((xr/cdelt1)+crpix1)
jj	=long((yr/cdelt2)+crpix2)
image(ii,jj)=255
image_sm0=smooth(image,3)
image_sm=smooth(image_sm0,5)
data(*,*,iday)=byte(image_sm)
tv,bytscl(data(*,*,iday),min=0,max=255)
wait,1
endfor

print,' '
print,'Simulated images and structures ____________ : ' 
index1	=index(0)
index2	=index(1)
data1	=data(*,*,0)
data2	=data(*,*,1)
coord_addindex,index1                  ;completes ephemerides R0,L0,B0,P in index
coord_addindex,index2

help,index1,index2,data1,data2
end