;+
; $Id: t_flimb.pro,v 1.1 2009/04/22 18:31:46 antunes Exp $
;
; Project   : developer area for STEREO SECCHI
;                   
; Name      : t_flimb
;               
; Purpose   : 
;               
; Explanation: 
;               
; Use       : IDL> 
;    
; Inputs    : 
;               
; Outputs   : 
;
; Keywords  : 
;
; Calls from LASCO : none 
;
; Common    : none 
;               
; Restrictions: none
;               
; Side effects: 
;               
; Category    : Pixon
;               
; Prev. Hist. : None.
;
; Written     : Sandy Antunes, NRL, Nov-Dec 2005
;               
; $Log: t_flimb.pro,v $
; Revision 1.1  2009/04/22 18:31:46  antunes
; first half of relocation.
;
; Revision 1.2  2009/04/21 19:08:51  antunes
; Updated up through 'r' with SSW headers.  Also minor John bug fixes.
;
; Revision 1.1  2009/03/31 14:57:01  antunes
; Re-org prepatory to checking into SSW.
;
; Revision 1.3  2006/11/16 18:28:06  antunes
; Added these back here, they were temporarily deleted.
;
; Revision 1.1  2006/04/11 15:55:06  antunes
; Massive re-org of cvs 'dev' preparatory to moving into solarsoft.
;
; Revision 1.2  2006/01/09 17:14:42  antunes
; More commented added.
; New 'super3drv' driver routine for getting user specs for reconstruction
; process now checked in.  Just type 'super3drv' to test.
;
;-            

pro t_flimb,ch,X,RR,d,Xo,fa,fb,fc,ImI0,an

; This subroutine returns the tangential and radial functions fa, fb and fc
; in the same format as the input coordinates X. These functions are
; to be used in the dsf and dsft functions for calculating electron density.
; For radii smaller than 1.1 solar radii, the functions are set to the value
; at 1.1 solar radii.
; 
; input   ch    - choice of limb darkening model
;         X     - coordinates of each vertex of the image (image units)
;         RR    - solar radius (image units)
;         d     - conversion factor
;         Xo  - image coordinates of sun center
; output  fa,fb,fc - extended sun functions Sigma_A, B, and C
;                    fa[i] is the value of fa at radius r (in image coordinates)
;                    where i=(r^2-rmin2)/(d*d)
;         ImI0     - the ratio Im/I0 where Im is mean solar intensity,
;                    and I0 is solar intensity at disc center.
;
;-----------------------------------------


NN=n_elements(X[0,*])

Ro  = double(RR)
Ro2 = Ro^2

if(d ne 1)then stop,'ERROR in flimb_07.pro' ; d is a fossil, I think

r2=d*((X[0,*]-Xo[0])^2+(X[1,*]-Xo[1])^2+(X[2,*]-Xo[2])^2)
r2 = reform(r2)
q=where(r2 le Ro,count)
if(count gt 0)then r2[q]=Ro2

case ch of

0: begin	; no limb darkening, point sun

   a0 = 1
   a1 = 0
   a2 = 0
   an = [a0,a1,a2]
   fa = Ro2/r2
   fb = Ro2/r2
   fc = 0
   ImI0=1.0
  end

1: begin	; polynomial limb darkening, point sun

   a1 =  0.93
   a2 = -0.23
   a0 = 1-a1-a2
   an = [a0,a1,a2]
   ImI0 = 1-a1/3-a2/2
   v  = ImI0*Ro2
   
   fa = v/r2
   fb = v/r2
   fc = 0

   end

2: begin	; polynomial limb darkening, extended sun

   a1 =  0.93
   a2 = -0.23
   a0 = 1-a1-a2
   an = [a0,a1,a2]
   ImI0 = 1-a1/3-a2/2

   s2 = Ro2/r2
   s  = sqrt(s2)
   c2 = 1-s2
   c  = sqrt(c2)
   c3 = c*c2
   q  = alog((1+s)/c)
   p  = (1-c)^2

   aa0 = (4-3*c-c3)/3
   aa1 = (5+s2-c2*(5-s2)*q/s)/8
   aa2 = 2*p*(c3+2*c2+8*c+4)/(15*s2)
   fa  = a0*aa0+a1*aa1+a2*aa2
     
   bb0 = c*s2
   bb1 = -(1-3*s2-c2*(1+3*s2)*q/s)/8
   bb2 = 2*p*(3*c3+6*c2+4*c+2)/(15*s2)
   fb  = a0*bb0+a1*bb1+a2*bb2

   fc = fa-fb

   end

else: stop,'ERROR flimb_07- ch = '+string(ch)
end

return
end