FUNCTION ellipse_maille , hdr, ihms
;+ ***********************************************************************
;NAME:
;	ELLIPSE_MAILLE
; PURPOSE:
;	Cette fonction calcule les parametres de l'ellipse section du lobe
;        a mi-puissance. Pour determiner cette ellipse on utilise la maille
;	 elementaire calculee dans maille_3.
;
; CATEGORY:
;	Traitement de fichiers NRH
; CALLING SEQUENCE:
;	
; INPUTS:
;       hdr     : header d'extension du fichier d'images FITS Nancay
;       
; COMMON BLOCKS:
;	RH	: communique avec les routines de lecture des donnees
;		  brutes RH	
; EXAMPLE:

;
; MODIFICATION HISTORY:

;		  lobe(3) : major minor axis in solar_radius and angle 
;			     in degrees of the major axis with Ox axis.
;
;	PSH, 2010/04/19: added "1" in call to INIT_MALAX_2D
;-**************************************************************************
;  Attention il faut d'abord avoir appele init_malax_2d qui initialise
;   le common malax.
; lecture dans le header des parametres necessaires a init_malax_2d
    date   = intarr(3)
    chdate = sxpar(hdr,'DATE-OBS')
    reads, chdate, date, format='(i4,1x,i2,1x,i2)'
    dat=date(0)
    date(0) = date(2)
    date(2) = dat    
    freq   = sxpar(hdr,'FREQ')
    chhmer = sxpar(hdr,'HMER')
    chdec  = sxpar(hdr,'DECL')
    ihmer  = intarr(4)
    reads, chhmer, ihmer, format='(i2,1x,i2,1x,i2,1x,i2)'
    idec   = intarr(4)
    reads,  chdec, idec,  format='(i3,1x,i3,1x,i3,1x,i3)'

    hmer      = fltarr(3)
    hmer(0:2) = float  (ihmer(0:2))
    hmer(2)   = hmer(2) + float(ihmer(3)) / 100.

    dec      = fltarr  (3)
    dec(0:2) = float   (idec(0:2))
    dec(2)   = dec (2) + float(idec(3)) / 100.

    sref = 'SOLEIL'		; 'soleil' ou 'autre'.
    rot_p = 1
    INIT_MALAX_2D,  date,  hmer,  dec,  sref, rot_p	; entree.
;   INIT_MALAX,  date,  hmer,  dec,  sref  ; meme resultat avec init_malax_2d 
;
    hmsc = mil_time(ihms)
    hmsc(3) = hmsc(3)/10.

    npi=128

    RH_MAILLE_3, npi,     freq,   hmsc, $           ; entree
              maille,  domega                       ; sortie
            ; entree: npi       nbre de pts sur l'image interferometrique.
            ;         frequence (MHz),
            ;         heure (h, m, s, c).
            ;         common MALAX, contenant entre autres new (defini par ie0,
            ;                        ie1, ie2) et nns (toujours 64).
;print,maille
;stop
    r1 = maille(2) - maille(0)		; unites  rayons solaires
    s1 = maille(3) - maille(1)
    r2 = maille(4) - maille(0)
    s2 = maille(5) - maille(1)

;   r1,s1 coordonnees heliog. du point de coordonnees interfero. (1,0)
;   r2,s2 coordonnees heliog. du point de coordonnees interfero. (0,1)
;print,r1,s1,r2,s2

; Taille du lobe :
;   En EW le lobe est decrit par un sinus cardinal qui s'annule au point 
;    d'abscisse 1 ; la 1/2 largeur a mi hauteur du lobe est  0.606
;   Pour le 2D le pave numerique donne en EW un lobe environ 3 fois plus large.
;
;   En NS, on surechantillonne un peu : on a 45 harmoniques sur 64 points
;                               
    coef_ew = 0.606 * 3.		; lobe 2D = lobe EW * 3
    coef_ns = 0.868		; = coef_ew*(64./45.)
    ae  = coef_ew		; Le rectangle circonscrit a l'ellipse a pour
    be  = coef_ns		;  1/2 cotes ae et be

; Dans le repere des coord. interferom. la section du lobe a mi-puissance est
;  une ellipse d 'axes Ox et Oy avec ae et be longueur des gd et petit axes..
; On cherche l'ellipse definie par ces 2 vecteurs : elle representera le lobe
; On passe du systeme interfero.(x,y) au systeme heliog.(X,Y) par la 
;  transformation :
;     X = r1x + r2y    
;     Y = s1x + s2y 
; La resolution du systeme donne :
;     x=(s2X-r2Y)/delta 
;     y=(r1Y-s1X)/delta
;     avec delta = r1s2-r2s1		; determinant du systeme

    delta = r1*s2-r2*s1

; Dans le repere des coord. interferom. (x,y)  l'ellipse  a pour equation:
;	x^2/ae^2 + y^2/be^2 = 1
;       ou ae et be sont les 1/2 axes de l'ellipse

; Dans le repere des coord. heliog. (X,Y)  l'ellipse a pour equation:	
;	AX2+BY2+CXY = 1
;  en remplacant x et y par leurs valeurs on trouve :
;	A = (s2^2 /ae^2 + s1^2/be^2) / delta^2
;	B = (r2^2 /ae^2 + r1^2/be^2) / delta^2
;	C = -2.*(r2*s2 /ae^2 + r1*s1/be^2) / delta^2

    A = (s2^2 /ae^2 + s1^2/be^2) / delta^2
    B = (r2^2 /ae^2 + r1^2/be^2) / delta^2
    C = -2.*(r2*s2/ae^2 + r1*s1/be^2) / delta^2
; determination des parametres de l'ellipse :
;    - g_a, p_a   : gd et petit axe de l'ellipse en rayons solaires
;    - teta       : angle en radians entre le gd axe et l'axe Ox.

; Pour trouver l'angle teta : On cherche le repere de coordonnees (XP,YP)
;  dans lequel l'equation de l'ellipse est :
;         AP*XP^2 + BP*YP^2 + CP*X*Y = 1 avec CP = 0  
;  cad gd axe et petit axe colineaires aux axes du repere. 
;   On effectue une rotation d'angle teta faisant passer du repere [ox,oy]
;   au repere OXP,OYP:
;              X = XP*costeta - YP*sinteta
;              Y = XP*sinteta + YP*costeta
; CP = 0  ==> On trouve teta = 0.5*atan ( C/(A-B))
;           et g_a = 1./(sqrt(A*costeta^2+B*sinteta^2+C*sinteta*costeta)
;	       p_a = 1./(sqrt(A*sinteta^2+B*costeta^2-C*sinteta*costeta)

    teta = 0.5*atan ( C/(A-B))
    costeta = cos(teta)
    sinteta = sin(teta)
    g_a = 1./(sqrt(A*costeta^2+B*sinteta^2+C*sinteta*costeta))
    p_a = 1./(sqrt(A*sinteta^2+B*costeta^2-C*sinteta*costeta))
;stop
    if (g_a lt p_a) then begin		; il faut echanger les axes
       jean_marc = g_a
	g_a = p_a
	p_a = jean_marc
	teta = teta + !pi / 2
	if (teta gt !pi) then  teta = teta - !pi
    endif
    teta = teta * !radeg 
;stop
; pour tracer tvellipse, g_a*solar_r, p_a*solar_r, 25, 25, teta 
    return,[g_a,p_a,teta]
    end