;+
;________________________________________________________________________________________________________________________
;
; p2sw_limb_profile : retourne un profil du limbe en polaires.
;________________________________________________________________________________________________________________________
;
; Appel : prof = p2sw_limb_profile(img, param, type, theta_size, rho_size, band, [gamma])
;
; Resultat : tableau de dimension [theta_size, 3] ou [theta_size, 4] (si le keyword gamma est present. prof[*, 0]
;            contient les valeurs de theta valides et prof[*, 1] contient les valeurs correspondantes de la distance au
;            centre du max du gradient. prof[*, 2] contient les points valides. prof[*, 3] contient les valeurs du gamma.
;
; Parametres : img : tableau 2x2. Image dans laquelle ont veut mesurer les coordonnees du centre et le rayon. 
;              param : tableau de trois elements contenant les coordonnes et le rayon suppose. 
;              type : chaine de caracteres decrivant le type de limbe (maximum ou minimum du gradient). Doit etre egal
;                     a 'min' ou 'max'.
;              theta_size : nombre d'echantillons le long d'une circonference. Une bonne valeur est 10000. 
;              rho_size : nombre d'echantillons le long d'un rayon. Une bonne valeur est 40.
;              band : tableau de deux elements contenant les bornes min et max (en rayons solaires) de la bande dans
;                     laquelle le limbe est suppose se situer.
;            
; Keywords : roll : angle antre le nord solaire et le haut de l'image. Permet le traitement des images tournees.
;            gamma : Tableau de deux elements contenant les valeurs min et max du contraste du limbe. si
;                    present, ne prend pas en compte dans le fit tous les points du limbe ne presentant pas un
;                    contraste compris entre les deux bornes.
;            holes : si present et different de 0, prend en compte les latitudes des trous coronaux pour le fit.
;                    0.0->0.1 CH | 0.1->0.4 | 0.4->0.6 CH | 0.6->0.9 | 0.9->1.0 CH.
;________________________________________________________________________________________________________________________
;
; Historique : 10/21/1998 : modification du critere de convergence.
;              28/01/2000 : ajout du keyword roll pour le traitement des images tournees.
;              23/09/2010 : call p2sw_polar_map with cubic to force
;                           cubic interpolation (dbs)
;________________________________________________________________________________________________________________________
;-
  
function p2sw_limb_profile, img, param, type, theta_size, rho_size, band, roll=roll, gamma=gamma, holes=holes

  holes_present = keyword_set(holes)
  gamma_present = keyword_set(gamma)

  min_r = param[2]*band[0]
  max_r = param[2]*band[1]

  x = 2.0*!pi*dindgen(theta_size)/theta_size

  a = p2sw_polar_map(double(img), theta_size, rho_size, param[0], param[1], -roll, -roll+360.0, min_r, max_r, missing=0, /cubic)

  derivee = (deriv(transpose(a)))[1:rho_size-2, *]
  a = a[*, 1:rho_size-2]
 
  prof = dblarr(theta_size) 
  if type eq 'min' then begin
    for k = 0l, theta_size-1 do begin
      dummy = min(derivee[*, k], m)
      prof[k] = m
    endfor
  endif else begin
    for k = 0l, theta_size-1 do begin
      dummy = max(derivee[*, k], m)
      prof[k] = m
    endfor
  end   
    
  dummy = rebin(double(a ne 0), theta_size, 1) eq 1.0
  if holes_present eq 0 then begin
    dummy[0:0.1*theta_size] = 0
    dummy[0.4*theta_size:0.6*theta_size] = 0
    dummy[0.9*theta_size:theta_size-1] = 0
  endif

  if gamma_present then begin
    mask = (replicate(1.0, theta_size) # dindgen(rho_size-2)) - rebin(prof, theta_size, rho_size-2, /sample)
    if type eq 'min' then mask = -temporary(mask)
    nb = (rho_size-2)*rebin(double((mask gt 0)), theta_size, 1)
    pgamma = ((rho_size-3-nb)/(nb>1)) * rebin(double((mask gt 0)*a), theta_size, 1) / rebin(double((mask lt 0)*a), theta_size, 1)
  endif

  if gamma_present eq 0 then $ 
    nonzero = dummy and (prof lt rho_size-3) and (prof gt 0) $
  else nonzero = dummy and (prof lt rho_size-3) and (prof gt 0) and (pgamma gt gamma[0]) and (pgamma lt gamma[1])
    
  prof = min_r + (prof + 1.5d) * (max_r - min_r) / rho_size
  if gamma_present then begin
    return, [[x], [prof], [nonzero], [pgamma]]
  endif else return, [[x], [prof], [nonzero]]

end