pro nlfff_alpha_2d,dir,file,vers,para
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : NLFFF_ALPHA
;
; Category    : Magnetic data modeling
;
; Explanation : Calculating nolinear force-free parameter ALPHA[NK,4]
;               averaged for every loop with two methods.
;
; Syntax      : IDL>magn_alpha,runfile,para
;
; Inputs      : runfile   = input file *_cube.sav, *_fit.sav, *_coeff.dat
;               para      = input parameters
;
; Outputs     ; runfile   = output file *_fit.sav
;			-->bzmap,x,y,field_lines,alpha,q_scale,angle,cpu,merit,iter
;
; History     : 16-May-2011, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

print,'_________________NLFFF_ALPHA__________________________'
nseg    =long(para[8])
hmax	=para[12]
ds	=para[13]

runfile=dir+file+vers
savefile=runfile+'_fit.sav'
restore,savefile      ;-->bzmap,x,y,field_lines,q_scale,angle,cpu,merit
help,dev_loop,dev_loop2
nx	=n_elements(x)
ny	=n_elements(y)
nz	=n_elements(z)
dim	=size(field_lines)
if (dim[0] le 2) then nf=1
if (dim[0] ge 3) then nf=dim[3]

coeff_file=runfile+'_coeff.dat'
readcol,coeff_file,m1,x1,y1,z1,a1
coeff   =[[m1],[x1],[y1],[z1],[a1]]

;_____________________FITTING POINT ARRAY_____________________________'
nh      =long(hmax/ds)
nfit    =nf*nh
x_fit   =fltarr(nfit,3)
v_fit   =fltarr(nfit,3)
for k=0,nf-1 do begin                           ;number of loops
 ind  =where(field_lines[*,0,k] ne 0,ns)	;2D coordinates (z=0)
 xl   =field_lines[0:ns-1,0,k]
 yl   =field_lines[0:ns-1,1,k]
 zl   =field_lines[0:ns-1,2,k]
 for ih=0,nh-1 do begin
  rh	=1.+hmax*float(ih)/float(nh-1)
  im    =long((ns-1)*0.5)
  ii    =k*nh+ih
  x_fit[ii,0]=xl[im]
  x_fit[ii,1]=yl[im]
  x_fit[ii,2]=sqrt(rh^2-xl[im]^2-yl[im]^2)
  if (ns gt 1) then begin
   v_fit[ii,0]=xl[im+1]-xl[im]
   v_fit[ii,1]=yl[im+1]-yl[im]
   v_fit[ii,2]=zl[im+1]-zl[im]
  endif
 endfor
endfor                                          ;for k=0,nf-1

;______________________CALCULATION OF ALPHA FIELDLINE____________________'
x0	=x_fit
x1	=x_fit
x2	=x_fit
y1	=x_fit
y2	=x_fit
z1	=x_fit
z2	=x_fit
x1[*,0]	=x1[*,0]-ds
x2[*,0]	=x2[*,0]+ds
y1[*,1]	=y1[*,1]-ds
y2[*,1]	=y2[*,1]+ds
z1[*,2]	=z1[*,2]-ds
z2[*,2]	=z2[*,2]+ds
nlfff_vector,coeff,x0,b_x0
nlfff_vector,coeff,x1,b_x1
nlfff_vector,coeff,x2,b_x2
nlfff_vector,coeff,y1,b_y1
nlfff_vector,coeff,y2,b_y2
nlfff_vector,coeff,z1,b_z1
nlfff_vector,coeff,z2,b_z2
dby_dx  =(b_x2[*,1]-b_x1[*,1])/(x2[*,0]-x1[*,0])
dbz_dx  =(b_x2[*,2]-b_x1[*,2])/(x2[*,0]-x1[*,0])
dbx_dy  =(b_y2[*,0]-b_y1[*,0])/(y2[*,1]-y1[*,1])
dbz_dy  =(b_y2[*,2]-b_y1[*,2])/(y2[*,1]-y1[*,1])
dbx_dz  =(b_z2[*,0]-b_z1[*,0])/(z2[*,2]-z1[*,2])
dby_dz  =(b_z2[*,1]-b_z1[*,1])/(z2[*,2]-z1[*,2])
jx	=(dbz_dy-dby_dz)
jy	=(dbx_dz-dbz_dx)
jz	=(dby_dx-dbx_dy)
bx	=b_x0[*,0]					
by	=b_x0[*,1]
bz	=b_x0[*,2]
wx	=bx^2						
wy	=by^2
wz	=bz^2
a       =(wx*(jx/bx)+wy*(jy/by)+wz*(jz/bz))/(wx+wy+wz)	;weighted by B^2

print,'______________ALPHA LOOP STATISTIC___________________________
alpha	=fltarr(nf,2)
for k=0,nf-1 do begin
 ind    =k*nseg+findgen(nseg)
 statistic,a[ind],a_avg,a_sig			;cube-calculated alpha
 alpha(k,0)=a_avg	
 alpha(k,1)=a_sig 	
endfor

;____________________CALCULATION OF ALPHA PHOTOSPHERE_______________'
amap	=fltarr(nx,ny)
bmap	=fltarr(nx,ny,4)
jmap	=fltarr(nx,ny)
for j=0,ny-1 do begin
 if (j mod 10) eq 0 then print,'j=',j
 xx	=x
 yy	=fltarr(nx)+y(j)
 zz	=sqrt(1.-xx^2-yy^2)			;curved photospheric surface
 x_fit	=[[xx],[yy],[zz]]
 x0	=x_fit
 x1	=x_fit
 x2	=x_fit
 y1	=x_fit
 y2	=x_fit
 z1	=x_fit
 z2	=x_fit
 x1[*,0]=x1[*,0]-ds
 x2[*,0]=x2[*,0]+ds
 y1[*,1]=y1[*,1]-ds
 y2[*,1]=y2[*,1]+ds
 z1[*,2]=z1[*,2]-ds
 z2[*,2]=z2[*,2]+ds
 nlfff_vector,coeff,x0,b_x0
 nlfff_vector,coeff,x1,b_x1
 nlfff_vector,coeff,x2,b_x2
 nlfff_vector,coeff,y1,b_y1
 nlfff_vector,coeff,y2,b_y2
 nlfff_vector,coeff,z1,b_z1
 nlfff_vector,coeff,z2,b_z2
 dby_dx  =(b_x2[*,1]-b_x1[*,1])/(x2[*,0]-x1[*,0])
 dbz_dx  =(b_x2[*,2]-b_x1[*,2])/(x2[*,0]-x1[*,0])
 dbx_dy  =(b_y2[*,0]-b_y1[*,0])/(y2[*,1]-y1[*,1])
 dbz_dy  =(b_y2[*,2]-b_y1[*,2])/(y2[*,1]-y1[*,1])
 dbx_dz  =(b_z2[*,0]-b_z1[*,0])/(z2[*,2]-z1[*,2])
 dby_dz  =(b_z2[*,1]-b_z1[*,1])/(z2[*,2]-z1[*,2])
 jx	=(dbz_dy-dby_dz)
 jy	=(dbx_dz-dbz_dx)
 jz	=(dby_dx-dbx_dy)
 bx	=b_x0[*,0]					
 by	=b_x0[*,1]
 bz	=b_x0[*,2]
 wx	=bx^2						
 wy	=by^2
 wz	=bz^2
 amap[*,j]=(wx*(jx/bx)+wy*(jy/by)+wz*(jz/bz))/(wx+wy+wz)	;weighted by B^2
 jmap[*,j]=jz
 bmap[*,j,0]=bx
 bmap[*,j,1]=by
 bmap[*,j,2]=bz
 bmap[*,j,3]=sqrt(bx^2+by+2+bz^2)
endfor

;_________________________LOOP STATISTIC_____________________________
alpha	=fltarr(nf,2)
for k=0,nf-1 do begin
 ind    =k*nseg+findgen(nseg)
 statistic,a[ind],a_avg,a_sig			;cube-calculated alpha
 alpha(k,0)=a_avg	
 alpha(k,1)=a_sig 	
endfor

;________________________TEST DISPLAY________________________________
test	=0
if (test eq 1) then begin
 for k=0,nf-1 do begin
  clearplot
  window,2,xsize=300,ysize=900
  !p.position=[0.1,0.7,0.95,0.95]
; plot,alpha__,yrange=amax*[-1,1]
  na	=n_elements(alpha__)
  oplot,findgen(na),fltarr(na)+a_avg,linestyle=2,thick=2
  oplot,alpha_,linestyle=2
  oplot,!x.crange,[0,0],linestyle=1
  xyouts_norm,0.5,0.9,string(a_avg,'(f5.2)')+'+'+string(a_sig,'(f5.2)'),char
  !noeras=1

  !p.position=[0.1,0.4,0.95,0.65]
  ind	=where(field_lines(*,2,k) ne 0,ns) 
  ineg	=where(alpha_ lt 0,ns) 
  plot,field_lines(ind,0,k),field_lines(ind,2,k),xrange=[-0.1,0.1],yrange=[0.9,1.1]
  if (ns ge 1) then oplot,field_lines(ineg,0,k),field_lines(ineg,2,k),thick=3

  !p.position=[0.1,0.1,0.95,0.35]
  plot,field_lines(ind,0,k),field_lines(ind,1,k),xrange=[-0.1,0.1],yrange=[-0.1,0.1]
  if (ns ge 1) then oplot,field_lines(ineg,0,k),field_lines(ineg,1,k),thick=3

  read,'continue?',yes
  if (yes eq 0) then stop
  wait,0.3
 endfor
endif

;_______________________SAVE ALPHA PARAMETER_________________________
save,filename=savefile,para,x,y,bzmap,field_loops,ns_loop,ns_field,$
      angle,dev_loop,iter,cpu,field_lines,z,bx,by,bz,b,a,merit,amap,$
      bmap,jmap,alpha,dev_loop2,amax
print,'parameters saved in file = ',savefile
end