pro fitscan,x,a,res,scl,wdt,sft,ofs
;PURPOSE Fits a Gaussian to a sequence of images (scan) at 
;	each pixel. Preferably rebin large images to smaller, lower
;	resolution arrays (unless you want to wait for hours...)
;INPUT	x: 1-d array(scanlength) describing scan sequence (e.g. 
;	    array of position angles in a polarizer scan).
;	a: 3-d array containing the images of the scan sequence
;	   (scan dimension first).
;OUTPUT	res: 3-d array of fitresult (same dimensions as a)
;	scl: 2-d "image" of resulting scale's
;	wdt: 2-d "image" of resulting e-folding widths
;	sft: 2-d "image" of resulting x-shifts
;	ofs: 2-d "image" of resulting offsets

common dimen,lnth,dlnth,area,darea,tarea,d
;Extract useful dimensions to be more efficient
  dim=size(a)
  lnth=dim(dim(0)+2) & dlnth=2*lnth
  if (dim(0) eq 1) then area=1 else area=dim(2)*dim(3)
  darea=2*area & tarea=3*area
  d=dim(1)
;Prepare "x,y,w" for curvefit (string out array of curves into one)
  xx=fltarr(lnth) & y=fltarr(lnth)
  xx(*)=x#replicate(1.,area)
  y(*)=a
  w=replicate(1.,lnth)
;Set initial guess for parameters - Edit to fit case <---- !!!!!!!!!
  p=fltarr(area*4)
  p(0:area-1)=replicate(1.,area)	; Scale
  p(area:darea-1)=replicate(2.5,area)	; Width (e-folding)
  p(darea:tarea-1)=replicate(5.,area)	; x-shift
;  p(tarea:*)=replicate(0.,area)	; Offset
;Now do the fit and get resulting array of curves...
  res=a
  res(*)=curvefit(xx,y,w,p)
;...and load resulting parameters.
  if (area eq 1) then scl=fltarr(1) else scl=fltarr(dim(2),dim(3))
  wdt=scl & sft=scl & ofs=scl
  scl(*)=p(0:area-1)
  wdt(*)=p(area:darea-1)
  sft(*)=p(darea:tarea-1)
  ofs(*)=p(tarea:*)
end

pro funct,x,p,f,pder
common dimen,lnth,dlnth,area,darea,tarea,d
  arg=fltarr(lnth) & f=fltarr(lnth) & pder=fltarr(lnth,area*4)
  sub=indgen(lnth)/d*long(lnth+d)+indgen(lnth) MOD d
  invwdt=replicate(1.,d)#p(area:darea-1)
  arg(*)=x-replicate(1.,d)#p(darea:tarea-1)
  arg=arg/invwdt
  invwdt=2./invwdt
  f=exp(-arg^2)
  pder(sub)=f
  f=f*(replicate(1.,d)#p(0:area-1))
  pder(sub+lnth*area*2)=f*invwdt*arg
  pder(sub+lnth*area)=pder(sub+lnth*area*2)*arg
  pder(sub+lnth*area*3)=1.
  f=f+replicate(1.,d)#p(tarea:*)
end