; Copyright (c) 1998-2008 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		AX_POL_FIT
;
;LAST CHANGED: ----------------------------------- 	19-Apr-08
;
; PURPOSE:
;	This function fits a curve of I(angle) to the expression for
; photoexcitation polarization dependence:
;        I = constant + amplitude(Angle-AngleRef)
; using CURVEFIT (IDL non-linear least squares fit)
; Equal weighting of all data points is used
;
; CATEGORY:
; Data analysis; called from axis2000 spectra~curefit~ax_pol_fit
; Used as part of fit to polarization maps
;
; CALLING SEQUENCE:
;	Result = AX_POL_FIT(CURV, MAX_ITER=MAX_ITER, PARAM=PARAM, SIGMA=SIGMA, TOL=TOL)
;
; INPUTS:
;  CURV:   (angle, intensity) 2-d array
;
; KEYWORDS:
; MAX_ITER:	if defined, the maximum number of iterations (default = 20)
; SIGMA: 	if defined, returns the standard definition
; SILENT:	if set, do not report results in axis or IDL log
; STATUS	if set, returns the status of the fit
;			    (0 = successful; 1 = divergent, 2 = non-convergent)
; PARAM;	if defined, gives & returns the (constant, amplitude, angle_reference) values
; TOL:    	if defined, the tolerance (default = 1e-8)
;
; OUTPUTS:
;  the fit function
;
; COMMON BLOCKS:
;	AXIS_COM	standard set of common blocks
;
; SIDE EFFECTS: none
;
; EXAMPLE: see IDL help file for example this was modelled from
;
; MODIFICATION HISTORY:
; (04-May-06 aph) first derived
; (25-Jan-07 aph) set default tolerance to 1e-6, max iter to 500
; (19-Feb-07 aph) change sign in fit function
; (19-Apr-08 aph) added /noderivative ot make ref angle fit better (from Martin Obst 25-Mar-08)
;-


PRO ang_funct,  angle, parm, pred, pder
; -------------------------------------------
; finction for non-linear curve fit to angle dependent linear dichroism
; parm(0) = constant
; parm(1)= amplitude
; parm(2)= angle_ref (in DEGREES)
;
angl = (3.1428/180.)*angle			; convert to radian
angl_ref = (3.1428/180.)*parm(2)		; convert to radian
del_ang = angl - angl_ref			; pre 19-feb-07  this was a SUM not a difference
    pred = parm(0) + parm(1)*cos(del_ang)*cos(del_ang)

;If the procedure is called with four parameters, calculate the
;partial derivatives.
  IF N_PARAMS() GE 4 THEN $
    pder = [[replicate(1.0, N_ELEMENTS(angle))],[cos(del_ang)*cos(del_ang)], $
              [-2.*parm(1)*cos(del_ang)*sin(del_ang)]]

END

FUNCTION AX_POL_FIT, curv, max_iter=max_iter, param=param, sigma=sigma, status = status, silent=silent, tol=tol

@axis_com
on_error,2

if not keyword_set(tol) then tol = 1e-6
if not keyword_set(max_iter) then max_iter=500

; ------ define default parameters

if not keyword_set(param) then param = [0.5,0.5,90]		; assumes application is spider silk

X=reform(curv(0,*))
Y=reform(curv(1,*))
; if wt is not defined (as here) no weighting is used
yfit = CURVEFIT(X, Y, wt, param, SIGMA, FUNCTION_NAME='ang_funct', $
        itmax=max_iter, iter=iter, status=status, tol=tol, /noderivative, /double)	;  added /noderivative

IF not keyword_set(silent) then begin
	CASE STATUS of
	0: axis_log, 'successful fit'
	1: axis_log, 'divergent fit'
	2: axis_log, 'non-convergent'
	ENDCASE
	axis_log, ' in ' + string(iter) + ' iterations'
	axis_log, 'Constant  = ' + string(param(0))
	axis_log, 'Amplitude = ' + string(param(1))
	axis_log, 'Ref_angle = ' + string(param(2))

ENDIF

return, yfit

END
