pro magn_mdidata,inputfile,para,fov
;+
; Project     : SDO, STEREO, SOHO
;
; Name        : MAGN_MDIDATA
;
; Category    : magnetic field moedling
;
; Explanation : reads MDI magnetogram data
;               and stores parameters in savefile
;
; Syntax      : IDL> magn_mdidata,filename,para
;
; Inputs      : inputfile= *._coeff.dat, *_cube.sav, *_field.sav,
;               para     = input parameters
;                          [ds,ngrid,nmag,nsm,qloop,nseg,acc,amax,hmax,dpix,
;                           niter,qzone,eps,qv,thresh,q_mag,halt,iopt]
;               inr      =1,2 (number of low Low&Lou data)
;
; Outputs     ; saves parameters in file *_cube.sav:  
;		--> x,y,z,bx,by,bz,b,a
;		saves parameters in file *_field.sav: 
;               --> bzmap,x,y,field_lines,alpha,q_scale,cpu,merit
;
; History     : 16-Feb-2012, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

print,'_______________MAGN_MDIDATA_________________________

;________________FIND LOOP FILE____________________________
search  =inputfile+'_loop.dat'
print,'Search for files with loop 3D coordinates : ',search
files	=findfile(search,count=nfiles)
if (nfiles eq 0) then stop,'No file found !'
loopfile=files(0)
readcol,loopfile,iloop,xl,yl,zl,sl
nf	=long(max(iloop)+1)
print,'Number of loop files   NFILES = ',nfiles
print,'Loop file            LOOPFILE = ',loopfile
print,'Number of loops            NF = ',nf 

;________________MAXIMUM LENGTH OF FIELDLINES______________
nsmax	=0
for i=0,nf-1 do begin
 ind	=where(iloop eq i,ns)
 nsmax	=nsmax>ns
endfor
print,'Max number of loop points  NS = ',nsmax

;________________READ MAGNETIC DATA________________________
search 	=inputfile+'_magn.fits'
print,'Search for files with magnetogram : ',search
magfile =file_search(search,count=nfiles)
if (nfiles eq 0) then stop,'No file found !'
image   =readfits(magfile(0),index_mag)
dateobs =sxpar(index_mag,'DATE_OBS')
nx      =sxpar(index_mag,'NAXIS1')
ny      =sxpar(index_mag,'NAXIS2')
crpix1  =sxpar(index_mag,'CRPIX1')
crpix2  =sxpar(index_mag,'CRPIX2')
cdelt1  =sxpar(index_mag,'CDELT1')
cdelt2  =sxpar(index_mag,'CDELT2')
pos0    =sxpar(index_mag,'SOLAR_P0')
rpix    =sxpar(index_mag,'R_SUN')               ;solar radius in pixels
dpix    =1./rpix
para[9]	=dpix	
print,'Number of magn files   NFILES = ',nfiles
print,'Loop file             MAGFILE = ',magfile(0)
print,'Date of observation   DATEOBS = ',dateobs
print,'Image dimensions        NX,NY = ',nx,ny
print,'Sun center position  CRPIX1,2 = ',crpix1,crpix2
print,'Pixel size           CDELT1,2 = ',cdelt1,cdelt2
print,'Solar radius (pixels)    RPIX = ',rpix
print,'Pixel size (solar radii) DPIX = ',dpix
print,'Rotation angle           POS0 = ',pos0


;________________EXTRACT MAGNETOGRAM_________________________
p0      =pos0*!pi/180.
x1	=fov[0]
x2	=fov[2]
y1	=fov[1]
y2	=fov[3]
i1	=long(crpix1+x1/dpix)
i2	=long(crpix1+x2/dpix)
j1	=long(crpix2+y1/dpix)
j2	=long(crpix2+y2/dpix)
ni	=i2-i1+1
nj	=j2-j1+1
x	=x1+dpix*findgen(ni)
y	=y1+dpix*findgen(nj)
bzmap	=image(i1:i2,j1:j2)
print,'FIELD OF VIEW             FOV = ',fov
print,'SUBIMAGE DIMENSION      NI,NJ = ',ni,nj

;________________ROTATION OF MDI IMAGE______________________
if (p0 ne 0.) then begin
 zmin	=min(bzmap)
 for j=0,nj-1 do begin
  y_     =fltarr(ny)+y(j)
  x_rot  = x*cos(p0)+y_*sin(p0)
  y_rot  =-x*sin(p0)+y_*cos(p0)
  i_     =long(rpix*x_rot+crpix1) 
  j_     =long(rpix*y_rot+crpix2) 
  stripe =image(i_,j_)
  ind0   =where(stripe eq zmin,nind0)
  if (nind0 ge 1) then stripe(ind0)=0.
  bzmap(*,j)=stripe
 endfor
endif

;________________DISPLAY______________________________________
window,0,xsize=ni,ysize=nj
loadct,0
!p.position=[0,0,1,1]
plot,[0,0],[0,0],xrange=[x1,x2],yrange=[y1,y2],xstyle=1,ystyle=1
tv,bytscl(bzmap)

;________________FIELD LINE ARRAY___________________________
loadct,5
field_lines=fltarr(nsmax,3,nf)
for i=0,nf-1 do begin
 ind	=where(iloop eq i,ns)
 field_lines(0:ns-1,0,i)=xl(ind)          ;x-coordinate loop
 field_lines(0:ns-1,1,i)=yl(ind)          ;y-coordinate loop
 field_lines(0:ns-1,2,i)=zl(ind)          ;z-coordinate loop
 oplot,xl(ind),yl(ind),color=128
endfor

;________________SAVE PARAMETERS______________________________
file1   =inputfile+'_cube.sav'
bz	=bzmap
bx	=fltarr(ni,nj)
by	=fltarr(ni,nj)
b	=bzmap
a	=fltarr(ni,nj)
nz	=1
z	=fltarr(nz)
save,filename=file1,x,y,z,bx,by,bz,b,a
print,'parameters saved in file = ',file1

file2	=inputfile+'_field.sav'
alpha	=fltarr(nf,2)
q_scale	=1.
cpu	=0.
merit	=[0.,0.,0.]
save,filename=file2,bzmap,x,y,field_lines,alpha,q_scale,cpu,merit
print,'parameters saved in file = ',file2
end