pro nlfff_image_stereo,files,adeg,lmin
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : NLFFF_FIELD
;
; Category    : Magnetic data modeling
;
; Explanation : Calculating NF magnetic field lines FIELD_LINES[NS,5,NF] 
;		magnetic field coefficients [B1,X1,Y1,Z1,A1] 
;
; Syntax      : IDL>nlfff_field,runfile,para
;
; Inputs      : runfile   = input file *_coeff.dat
;               para      = input parameters
;
; Outputs     ; runfile   = output files *_fit.sav
;
; History     : 16-May-2011, Version 1 written by Markus J. Aschwanden
;               30-Jan-2012, Version 2 released to SSW
;               29-JAN-2015, im --> np/2 (footpoint --> midpoint)
;
; Contact     : aschwanden@lmsal.com
;-

;______________________READ DATA_____________________________________
restore,files(5)	;-->para,bzmap,xaxis,yaxis,ns_field,field
nsmax	=long(para[7])
dfoot	=para[2]
dsmag	=para[3]
nseg	=para[8]
hmax	=para[12]
ds      =para[13]
xf1	=para[17]
yf1	=para[18]
xf2	=para[19]
yf2	=para[20]
asep	=!pi*(adeg/180.)		;half of AB separation
nx      =n_elements(xaxis)
ny      =n_elements(yaxis)
dpix    =float(xf2-xf1)/float(nx-1)
dim	=size(field)
if (dim(0) eq 2) then nf=1
if (dim(0) eq 3) then nf=dim(3)
xhalf	=(xf2-xf1)/2.
xmid	=(xf1+xf2)/2.
ymid	=(yf1+yf2)/2.
lat	=asin(ymid)
rho	=cos(lat)
l1	=asin(xmid/rho)			;longiutde of AR in STEREO-A
l2	=asep-l1			;longitude of AR in STEREO-B
xmid2	=-rho*sin(l2)			;xcoordinate of AR in STEREO-B
image1	=fltarr(nx,ny)
image2	=fltarr(nx,ny)
xaxis1	=xaxis
yaxis1	=yaxis
xaxis2  =xmid2-xhalf+(xaxis1-xaxis1(0))    ;x-axis of AR in STEREO-B
yaxis2	=yaxis
field2  =fltarr(nsmax,3,nf)

;______________________MAGNETIC FIELD LINES______________________
window,2,xsize=1024,ysize=1024
!p.position=[0.05,0.05,0.95,0.95]
!x.range=minmax(xaxis)
!y.range=minmax(yaxis)
!x.style=1
!y.style=1
plot,[1,1],[1,1]
for k=0,nf-1 do begin
 image3	=fltarr(nx,ny)
 image4	=fltarr(nx,ny)
 ns	=ns_field(k)
 xe	=field(0:ns-1,0,k)
 ye	=field(0:ns-1,1,k)
 oplot,xe,ye
;______________________STEREOSCOPIC ROTATION____________________
 ze	=field(0:ns-1,2,k)
 r	=sqrt(xe^2+ye^2+ze^2)
 h	=r-1.
 sin_a	=ye/(1.+h)
 cos_a	=sqrt(1.-sin_a^2)
 rho	=cos_a*(1.+h)
 lon_e	=asin(xe/rho)			;longitude of AR in Earth view 
 lon_a	=asep-lon_e			;longitude of AR in STEREO  
 xa	=-rho*sin(lon_a)		;xcoordinate of AR in STEREO
 ya	=ye
 za	=sqrt(r^2-xa^2-ya^2)
 field2(0:ns-1,0,k)=xa
 field2(0:ns-1,1,k)=ya
 field2(0:ns-1,2,k)=za
;_______________________CREATE EUV IMAGE________________________
 iy	=(ye-yaxis1(0))/dpix 
 ixe	=(xe-xaxis1(0))/dpix 
 ixa	=(xa-xaxis2(0))/dpix
 temp	=2.0	;MK
 lambda =(temp*50./696.)/2.
 image3(ixe,iy)=exp(-h/lambda)
 image4(ixa,iy)=exp(-h/lambda)

 w0	=(0.0015/dsmag)		;standard width 3 HMI pixels 0.0005*3
 w	=w0*alog10(sqrt(ns))	;Gaussian half width
 nk	=long(w*6+1)		;extent of wings
 nk2	=nk/2			;center position
 kernel	=fltarr(nk,nk)
 for i=0,nk-1 do begin
  for j=0,nk-1 do begin
   r	=sqrt((i-nk2)^2+(j-nk2)^2)
   kernel(i,j)=exp(-r^2/(2.*w^2))
  endfor
 endfor
 image1=image1 + convol(image3,kernel)
 image2=image2 + convol(image4,kernel)
endfor

;________________________DISPLAY EUV IMAGE_____________________
window,2,xsize=2.*nx+1,ysize=ny
statistic,image1,avg,sig
z1	=0
z2	=sig*3
image2(0,*)=255
tv,bytscl(image1,min=z1,max=z2),0,0
tv,bytscl(image2,min=z1,max=z2),nx+1,0

;_______________________SAVE MODEL_____________________________
save,filename=files(0),para,bzmap,xaxis,yaxis,ns_field,field,field2,$
	image1,image2,xaxis1,yaxis1,xaxis2,yaxis2
print,'parameters saved in file = ',files(0)
end