PRO RH_ELLIPSE_MAILLE, ellipse, freq, hmsc, maille = ma ;+ *********************************************************************** ;NAME: ; RH_ELLIPSE_MAILLE() ; PURPOSE: ; Cette fonction calcule les parametres de l'ellipse qui approche la ; section du lobe a mi-puissance. Pour determiner cette ellipse on ; utilise la maille elementaire calculee dans RH_MAILLE_3. ; La maille elementaire est calculee dans cette fonction ; ou bien est fournie en entree par le KEYWORD MAILLE.(dans le cas ; ou elle est calculee en amont dans le programme appelant) ; Attention il faut d'abord avoir appele INIT_MALAX qui initialise ; le common malax. ; ; CATEGORY: ; Traitement de fichiers NRH ; CALLING SEQUENCE: ; OUTPUTS: ; ELLIPSE tableau contenant : ; major axis, minor axis (in solar_radii), and angle in ; degrees of the major axis with Ox axis. ; ; OPTIONAL INPUTS: ; (si non founis, le keyword MAILLE doit etre present.) ; FREQ frequence (MHz) ; HMSC heure (heures, minutes, secondes, centiemes) ; ; OPTIONAL KEYWORD INPUT: ; MAILLE tableau contenant la maille elementaire quand il a ete ; calcule en amont dans le programme appelant. ; COMMON BLOCKS: ; RH : communique avec les routines de lecture des donnees ; brutes RH ; EXAMPLE: ; ; MODIFICATION HISTORY: ;-************************************************************************** ; DEBUT ESPACE COMMENTAIRES ; Rappel : si on appelle ; . r1, s1 coordonnees heliog. du point de coordonnees interfero. (1,0) ; . r2, s2 ------------------------------------------------------ (0,1) ; . maille le tableau de sortie de RH_MAILLE_3 donnant les coordonnees helio- ; graphiques (en rayons solaires) des 3 points definissant la maille ; elementaire, ; on a : ; r1 = maille(2) - maille(0) ; s1 = maille(3) - maille(1) ; r2 = maille(4) - maille(0) ; s2 = maille(5) - maille(1) ; Definition d'une ellipse approchant la section du lobe a mi_puissance ; Le lobe 2D a une forme compliquee du fait des "extensions filiformes" du ; "pave central". Une bonne methode consisterait a ajuster par moindres ; carres une gaussienne elliptique au lobe reel dans l'espace heliographi- ; que et a calculer les axes de sa section a mi-puissance. Une telle gaus- ; sienne elliptique serait representative du lobe "propre" obtenu par ; clean. ; On recule par flemme devant cet effort et on se contente de la demarche ; approchee suivante : ; En 1D : ; . en EW le lobe en coordonnes interferometriques est decrit par un ; sinus cardinal qui s'annule au point d'abscisse 1. Sa 1/2 largeur ; a mi-hauteur du lobe est 0.606 ; . en NS, on surechantillonne un peu : 128 points est le minimum de ; Shannon pour les 64 harmoniques EW et il n'y en a que 45 en NS. ; La largeur a mi-hauteur du lobe NS en coordonnees interferometri- ; ques est donc 0.606 * 64 / 45 = 0.862 ; En 2D le pave central seul donne un lobe de 1/2 largeurs a mi-hauteur : ; . 0.606 * 64 / 16 = 2.42 en EW ; . 0.862 * 45 / 23 = 1.69 en NS. ; Il s'y superpose des "tumnnels" 1D EW et NS plus etroits provenant ; des "extensions filiformes" correspondant aux antennes d'extension et ; ayant des 1/2 largeurs a mi-hauteur egales a celles des lobes 1D. ; On admet que la gaussienne elliptique ajustee au lobe par moindres carres ; a une section a mi-hauteur dont les 1/2 axes sont intermediaires entre ; les dimensions correspondant au pave central seul et celles correspon- ; dant aux cas 1D : ; . EW : entre 0.606 et 2.42 ; . NS : entre 0.862 et 1.69 ; La recette de cuisisne adoptee pour calculer ces valeurs intermediaires ; est de faire un moyenne ponderee par les nombres d'harmoniques coores- ; pondants: 23 * 16 = 368 pour le pave central seul, 48 pour les harmo- ; niques EW avec E0, E1, E2, 45+16=61 pour les harmoniques NS purs et ; ceux obtenus avec NS45 et les EW. ; On adopte donc les dimensions suivantes en EW et NS dans l'espace des ; coordonnees interferometriques sont donc : ; . (0.606 * 48 + 2.42 * 368) / (48 + 368) = 2.21 en EW ; . (0.862 * 61 + 1.69 * 368) / (61 + 368) = 1.57 en NS ; ; Calcul des axes de l'ellipse transformee dans l'espace des coordonnees helio- ; graphiques ; On passe du systeme interfero. (x, y) au systeme heliog. (X, Y) par la ; transformation : ; X = r1 * x + r2 * y (1) ; Y = s1 * x + s2 * y ; ou x et y sont en "canaux" et X et Y en rayons solaires. ; La resolution du systeme (1) donne : ; x = (s2 * X - r2 * Y) / delta (2) ; y = (r1 * Y - s1 * X) / delta ; avec delta = r1 * s2 - r2 * s1 (determinant du systeme). ; ; Equations de l'ellipse dans les espaces interferometriques et heliographiques ; . dans le repere interferometrique (x,y), dont les axes coincident avec ; ceux de l'ellipse, l'ellipse a pour equation: ; x^2 / ae^2 + y^2 / be^2 = 1 (3) ; ou ae et be sont les 1/2 axes de l'ellipse ; . dans le repere heliographique (X,Y) l'ellipse est a-priori inclinee ; sur les axes et a la forma plus generale : ; A * X^2 + B * Y^2 + C * X * Y = 1 (4) ; Les coefficients A, B, C s'obtiennent en substituant les valeurs de x et ; y en fonction de X et Y donnes par (2) dans l''equation (3). On trouve: ; A = (s2^2 / ae^2 + s1^2 / be^2) / delta^2 (5) ; B = (r2^2 / ae^2 + r1^2 / be^2) / delta^2 ; C = -2 * (r2 * s2 / ae^2 + r1 * s1 / be^2) / delta^2 ; ; Calcul des parametres de l'ellipse dans l'espace heliographique ; . g_a, p_a : grand axe et petit axe en rayons solaires, ; . teta : angle en radians entre le gd axe et l'axe OX. ; Pour trouver l'angle teta on cherche le repere principal (XP, YP) des coor- ; donnees, dans lequel l'equation de l'ellipse est : ; AP * XP^2 + BP * YP^2 + CP * X * Y = 1 avec CP = 0 . ; On effectue une rotation d'angle teta faisant passer du repere (OX, OY) ; au repere (OXP, OYP) : ; X = XP * costeta - YP * sinteta (6) ; Y = XP * sinteta + YP * costeta ; En subsituuant les valeurs de X et Y donnees par (6) dans l'equation (4) ; et en imposant la condition CP = 0, on trouve : ; teta = 1/2 * atan( C / (A - B) ) ; 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) ; FIN ESPACE COMMENTAIRES ;______________________________________________________________________________ ; Debut de la partie active if ( not keyword_set(ma)) then begin NPI = 128 ; nbre de points adoptes sur l'image interferometrique ; (ne pas changer) 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). Il a ete prealable- ; ment rempli par un appel a INIT_MALAX. ; sortie : maille est en rayons solaires ma = maille endif maille = ma r1 = maille(2) - maille(0) s1 = maille(3) - maille(1) r2 = maille(4) - maille(0) s2 = maille(5) - maille(1) ; print,r1,s1,r2,s2 ; stop ; Definition des axes de l'ellipse dans l'espace interferometrique ae = 2.21 ; en canaux interferometriques. be = 1.57 ; ---------------------------- delta = r1 * s2 - r2 * s1 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 ;print,'A B C ',A,B,C 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 jean_marc = 0 endif teta = teta * !radeg ;stop ; pour tracer sous IDL faire ; ; TVELLIPSE, g_a * solar_r, p_a * solar_r, 25, 25, teta ; ou solar_r est le nbre de pixels par rayon solaire (voir le ; mot-cle SOLAR_R des fichiers images Fits nrh2******.fts) ellipse = [g_a, p_a, teta] return end ; fin de ELLIPSE_MAILLE.