;---------------------------------------------------------------------
;
;                      ALLEN_FUN_O6.PRO
;
;---------------------------------------------------------------------
;
; Purpose
;
; This function computes the center-to-limb variation
;
; Called by FOR_UV_STOKES
;
; CALLING SEQUENCE:
;   allen_fun(mu)


function allen_fun_o6, mu

; mu = cos(zeta)= sqrt(1.0 - rho^2) for far-away placed observer (Data from Withbroe 1970)
;x = [1.0 , 0.99932963, 0.99731581, 0.99378125, 0.9878233, 0.98195013, 0.97383179, 0.96538475, 0.95330476, 0.94037119, 0.92443898, 0.91056769, 0.89381294, 0.87448642, 0.84844722, 0.82552196, 0.79927158, 0.7716396 , 0.74133063, 0.7095054 , 0.67618704, 0.64318597, 0.6050072 , 0.5535439 , 0.49938532, 0.44838606, 0.39602437, 0.34644578, 0.30671935, 0.26590224, 0.21777052, 0.17255434, 0.14453287, 0.0947615]

; o6_clv ratio from Withbroe 1970 Fig.5
;y = [1.0, 1.024, 1., 0.976, 1.048, 1.096, 1.096, 1.096, 1.096, 1.168, 1.168, 1.168, 1.24, 1.312, 1.36, 1.384, 1.456, 1.552, 1.6, 1.72, 1.816, 1.936, 2.032, 2.224, 2.464, 2.632, 2.896, 3.208, 3.568, 4., 3.69388, 3.38776, 3.08163, 2.83673]


; Limb darkening equation from Allen's book is a 2-deg polynomial which has been used for FORCOMP/uv.f in case of visible spectropolarimetry
;I_by_I0 = 1. - u2 - v2 + (u2 * mu) + (v2 * mu^2) where mu=cos(theta)


; Limb darkening for white-light uses a simpler equation from Ch. 6 eq. 13 (Billings 1966) which is same as Allen if v2=0 which has been used in OBSERVABLES/for_cfun.pro and for_bfun.pro

; The fit below is specific to O VI and requires a higher-order polynomial fit
; The variable ulimb will be used to scale this fit for the purposes of for e.g studying solar cycle effects and in general as a diagnostic tool for understanding the potential impact of symmetry breakigng

; The default choice of ulimb = 0 will turn off the center-to-limb variation effect.

;coeff = poly_fit(x, y, 11, /double)

c0 = -2.526
c1 = 49.94
c2 = -154.89
c3 = 140.86
c4 = 67.29
c5 = -174.31
c6 = 73.04
c7 = 1.035
c8 = 0.447
;c9 = -2379
;c10 = -5.357e04

;I_by_I0 = coeff[0] + (coeff[1] * mu) + (coeff[2] * mu^2) + (coeff[3] * mu^3) + (coeff[4] * mu^4) + (coeff[5] * mu^5) + (coeff[6] * mu^6) + (coeff[7] * mu^7) + (coeff[8] * mu^8) + (coeff[9] * mu^9) + (coeff[10] * mu^10) + (coeff[11] * mu^11) 

;I_by_I0 = c0 + (c1 * mu) + (c2 * mu^2) + (c3 * mu^3) + (c4 * mu^4) + (c5 * mu^5) + (c6 * mu^6) + (c7 * mu^7) + (c8 * mu^8) + (c9 * mu^9) ;+ (coeff[10] * mu^10) + (coeff[11] * mu^11) 

I_by_I0 = c0 + (c1 * mu) + (c2 * mu^2) + (c3 * mu^3) + (c4 * mu^4) + (c5 * mu^5) + (c6 * mu^6) + c7 * exp(c8 * mu)

;print, "11th deg polynomial", I_by_I0
return, I_by_I0

end