function get_noise, data1d

a=convol(data1d, [-1.,2.,-1.])
n=n_elements(a)
return, stdev(a[1:n-2])/sqrt(6.)
end

function niris_cogmag, wv, idata, vdata, wlc,  wv1=wv1, wv2=wv2
;
;+
;   NAME:  ch_cogmag
;   PURPOSE:  determine the line-of-sight field from Stokes I and V profiles
;          using the cener of gravity method
;   Calling Sequence:
;
;          b  = ch_cogmag(wv, idata, vdata, wlc, line=line)
;   INPUTS:
;           wv        wavelength (1D) in angstrom
;           idata, vdata        Stokes I and V  data (2D) normalized by the continuum intensity
;                      idata(*,0) --> specral profiles at 0-th point
;                      idata(0,*) ---> spatial profiles at 0-th wavelength
;   Optional output
;          wlc     center-of-gravity of wavelength
;  OUTPUTS:
;          b  an array of line-of-sight field
;  HISTORY
;         2007 June :  first coded by J Chae
;         2007 July :  dispersion "dwl" is replaced by wavelength array "wl"
;                      wavelength center is obtained as an optional output
;-
;common _MEsinglet,  geff, lambda_rest ; Lande g-factor and rest wavelength in angstrom

geff=3.0
lambda_rest=15648.5 ;

nwv = n_elements(idata[*,0])
npoint =n_elements(idata[0,*])


wv1=fltarr(npoint) & wv2=fltarr(npoint)

for k=0, npoint-1 do begin
if max(abs(vdata(*,k))) gt 0.03 then begin  
  
weight1=-1.*(vdata(*,k)<0);(1-(idata(*, k)+vdata(*,k)));was 1-(idata(*, k)+vdata(*,k))
weight2=vdata(*,k)>0;(1-(idata(*,k)-vdata(*,k)))
endif else begin
 weight1= (1-(idata(*, k)+vdata(*,k)))
 weight2= (1-(idata(*, k)-vdata(*,k)))
  
endelse

cutlevel1=  (abs(min(weight1)))
cutlevel2=  (abs(min(weight2)))
wv2[k] = total(weight1*wv)/total(weight1)
wv1[k] = total(weight2*wv)/total(weight2)
endfor

coeff=4.67e-13*lambda_rest^2*geff

B = (wv2-wv1)/(2*coeff);2. special coeff.

wlc = 0.5*(wv2+wv1)


return, B
end

pro niris_MEsinglet, x, par, f

;common ch_MEsinglet,  geff, lambda_rest ; Lande g-factor and rest wavelength in angstrom
geff=3.0
lambda_rest=15648.5 ;

B       = par[0] ; magnetic field strength in Gauss
theta   = par[1] ; inclination in radian
chi     = par[2] ;  azimuth in radian
eta0    = par[3] ;  ratio of line center opacity to the continuum one, fixed to 1.3
dlambdaD= par[4] ; Doppler width in angstrom, fixed to 0.15
a       = par[5] ; damping parameter, fixed to 1
lambda0 = par[6] ; line center in A
B0      = par[7] ; zeroth order of the source function
B1      = par[8] ; first order of the source function

lambda = x[0:n_elements(x)/4-1]
v= (lambda-lambda0)/dlambdaD

vb = geff*(4.67e-13*lambda_rest^2*B)/dlambdaD  ;

ch_voigt,  a, v+vb, phib, psib ;, phibda, phibdu, psibda, psibdu
ch_voigt,  a, v,    phip, psip ;, phipda, phipdu, psipda, psipdu
ch_voigt,  a, v-vb, phir, psir ;, phirda, phirdu, psirda, psirdu

factor=1./sqrt(!pi)
phib=phib*factor  & psib=psib*factor
phip=phip*factor  & psip=psip*factor
phir=phir*factor  & psir=psir*factor

st=sin(theta)
st2 = st^2
ct = cos(theta)



etaI =  1 + 0.5*eta0*(phip*st2+0.5*(phib+phir)*(1+ct^2))
etaQ = eta0*0.5*(phip-0.5*(phib+phir))*st2*cos(2*chi)
etaU = eta0*0.5*(phip-0.5*(phib+phir))*st2*sin(2*chi)
etaV = eta0*0.5*(phir - phib)*ct

rhoQ = eta0*0.5*(psip -0.5*(psib+psir))*st2*cos(2*chi)
rhoU = eta0*0.5*(psip -0.5*(psib+psir))*st2*sin(2*chi)
rhoV = eta0*0.5*(psir - psib)*ct

Delta = etaI^2*(etaI^2-etaQ^2-etaU^2-etaV^2+rhoQ^2+rhoU^2+rhoV^2) $
  - (etaQ*rhoQ+etaU*rhoU+etaV*rhoV)^2

I = B0 + B1*etaI*(etaI^2+rhoQ^2+rhoU^2+rhoV^2)/Delta
Q = -B1*(etaI^2*etaQ+etaI*(etaV*rhoU-etaU*rhoV) $
      +rhoQ*(etaQ*rhoQ+etaU*rhoU+etaV*rhoV))/Delta
U = -B1*(etaI^2*etaU+etaI*(etaQ*rhoV-etaV*rhoQ)  $
      +rhoU*(etaQ*rhoQ+etaU*rhoU+etaV*rhoV))/Delta
V = -B1*(etaI^2*etaV+etaI*(etaU*rhoQ-etaQ*rhoU)   $
      +rhoV*(etaQ*rhoQ+etaU*rhoU+etaV*rhoV))/Delta

f=[I,  Q, U, V]


end


function niris_mefit, x, data, par, chisq=chisq, display=display
n=n_elements(x)/4

weight=replicate(get_noise(data[0:n-1]), n)
;stop
for k=1, 3 do weight=[weight, replicate(get_noise(data[n*k:n*k+n-1]), n)]
weight=1./weight^2


;stop
 parinfo = replicate({value:0., fixed:0, limited:[0,0], $
                         limits:[0.,0]}, 9)
 ; parinfo[fixed].fixed = 1
  parinfo[[3,4,5,7,8]].limited[0] = 1
  parinfo[[3,4,5,7,8]].limits[0]  = 0.;  [0., 0, 0.1, 0.1,

;TEMP
  parinfo[[3, 5]].limited[1]=[1,1] & parinfo[[3,5]].limits[1]=[20, 5];was, [1,1], [10,5]
  parinfo[0].limited=[1,1]
  parinfo[0].limits=[0,6000]
  parinfo[1].limited=[1,1]
  parinfo[1].limits=[0, !pi]
  parinfo[2].limited=[1,1]
  parinfo[2].limits=[0, !pi]
  ;parinfo[3].fixed=1
  parinfo[4].fixed=1
  ;parinfo[5].fixed=1

  ;parinfo[3].limited[1]=1 & parinfo[3].limits[1]=5.

  parinfo[*].value = par ;  [5.7D, 2.2, 500., 1.5, 2000.]

if display then begin
;print, par
niris_MEsinglet, x, par, f
endif

result=mpcurvefit(x, data, weight, par, sigma, function_name='niris_MEsinglet' ,  quiet=1, $
chisq=chisq, noderiv=1, itmax=100, tol=1.e-4, parinfo=parinfo)

if keyword_set(display) then begin
window, 2, xs=800, ys=800
if n_elements(result) gt 1 then begin
for k=0, 3 do begin
if k eq 0 then yr=[0.5, 1.2] else yr=[-0.1, 0.15]
plot, x[0:n-1], data[n*k:n*(k+1)-1],  $
psym=1, xst=1, yst=1, yr=yr, ytitle=(['I', 'Q', 'U', 'V']+'/I!d0!n')[k], xtitle='Wavelength (A)'
oplot,x[0:n-1], f[n*k:n*(k+1)-1], linest=2, thick=2
oplot,x[0:n-1], result[n*k:n*(k+1)-1], linest=0, thick=2

endfor
endif

;print, par, chisq/(n*4)
;print, 'blos=', par[0]*cos(par[1]), ',  bt=', par[0]*sin(par[1])

endif
return, result
end

pro init_par, x, data, par, display=display

dwl_line=1.5;was, 1.5
dwl_core=0.15
btcoeff=5000.
n=n_elements(x)/4

continuum =median((data[0:n-1])[where(abs(x[0:n-1]) ge dwl_line )] )

line=where(abs(x[0:n-1]) le dwl_line)
blos=niris_cogmag((x[0:n-1])[line], (data[0:n-1])[line]/continuum, $
     (data[3*n:4*n-1])[line]/continuum)


core=where(abs(x[0:n-1]) le dwl_core ) ; or abs(x[0:n-1]-wlc[0]+0.4) le 0.1 )
icore=((data[0:n-1])[core])
qicore=((data[n:2*n-1])[core]) /icore
uicore=((data[2*n: 3*n-1])[core])/icore
bt=btcoeff*sqrt(sqrt(mean(qicore^2+uicore^2)))


;if keyword_set(display) then  print, 'blos=', blos[0], 'bt=', bt[0]

Bfield=sqrt(blos^2+bt^2) ; magnetic field strength in Gauss


theta= atan(bt, blos) ; inclination in radian

chi=atan(-uicore, -qicore)
chi=0.5*median(chi+2*!pi*(chi lt 0.))


eta0    = 1.3 ;  ratio of line center opacity to the continuum one
dlambdaD= 0.15 ; Doppler width in angstrom
adamp   = 1. ; damping parameter
lambda0 = 0. ; line center in A
 ; zeroth order of the source function
 ;/////////////////
B1 = 4*(continuum-min(icore)) ; first order of source function
B0=continuum-B1

par=[Bfield, theta, chi, eta0, dlambdaD,adamp, lambda0, B0, B1]
;       0      1     2     3      4       5       6      7   8

end


;MAIN  _______________________________________

function main_test, a, h

display=0

platform='win'
  set_plot, platform

sza=size(a)
nx=sza(1) & ny=sza(2) & nw=sza(3)
print, sza
b=a

lambda_rest=fxpar(h, 'refwv')
dlambda=findgen(fxpar(h, 'NAXIS3'))*fxpar(h, 'INCWV')+fxpar(h, 'STARTWV')
!p.multi=[0, 2,2]
s=where(abs(dlambda) le 2.5, n)
factor=median(b[*,*,40<(nw-1),0])

if 1 then begin
bmag=fltarr(nx, ny) & inclination=bmag & azimuth=bmag
chisqa=fltarr(nx, ny)
inta=fltarr(nx,ny)
parmap=fltarr(nx,ny,9)
fit_array=fltarr(nx,ny,n_elements(s),4)
endif

t1=systime(/s)
for xpos=0, nx-1 do for ypos=0, ny-1 do begin
;  for xpos=484, 484 do for ypos=298, 298 do begin
  ;for xpos=310, 310 do for ypos=324, 324 do begin

x=[dlambda[s], dlambda[s], dlambda[s], dlambda[s]]
b1=b[xpos:xpos,ypos:ypos, s, *]
cont=where(abs(x[0:n-1]) ge 1.5)

data=reform(b1)/factor

; adjustment of non-zero bias in Q, U and V?????
for k=1, 3 do data[k*n:(k+1)*n-1]=data[k*n:(k+1)*n-1]-median((data[k*n:(k+1)*n-1])[cont])

init_par, x, data, par, display=display
result=niris_mefit( x, data,  par, chisq=chisq, display=display)
if (n_elements(result)-1) gt 0 then fit_array[xpos, ypos, *,*]=reform(result, n_elements(s), 4)



if n_elements(par) eq 9 then begin
bmag[xpos,ypos]=par[0]
inclination[xpos, ypos]=par[1]
azimuth[xpos,ypos]=par[2]
chisqa[xpos, ypos]=chisq/(n*4)
inta[xpos, ypos]=par[7]+par[8]
parmap[xpos, ypos, *]= par
endif
if ypos eq 0 and xpos mod 50 eq 0 then print, 'Processed xpos of '+string(xpos)
endfor


t2=systime(/s)
print, 'Computing time=', t2-t1, ' s'

;window, xs=nx*2, ys=ny*2
;loadct, 13
;device, decomposed=0
;
;tv, bytscl(parmap[*,*,1], 0, !pi), 0
;tv, bytscl(parmap[*,*,2], 0, !pi), 1
;loadct, 3
;tv, bytscl(parmap[*,*,0], 0, 3000), 2
;loadct, 6
;tv, bytscl(-parmap[*,*,6], -0.3, 0.3), 3

loadct, 0
stop
return, parmap
end