; Liste des procedures : ; SIMUL_VISIB_1, ; SIMUL_VISIB_2, ; SIMUL_HARM_N, ; SIMUL_RH, qui produit le lobe theorique et l'image d'un soleil a partir ; d'une repartition d'harmoniques comportant les antennes anti- ; aliasing. ; RH_SIMULATION, bidon, apple'e pour compiler le tout. ;______________________________________________________________________________ PRO T_SIMUL_VISIB_1 ; Creation 9 jan 2006, pour essai de lancement de SIMUL_VISIB_1 direct_auto = 1 i_champ_uv = 1 SIMUL_VISIB_1, $ ; dans fichier RH_SIMULATION.PRO ; entree direct_auto, $ i_champ_uv , $ ; bidon ; sortie objet_i, $ ; Cygne centre et defini sur un champ ~ 22' d'arc dans ; SIMUL_VISIB_1 (fichier CALIB_SUB) par 512 points ; pour etre bien decrit : ; . Cygne dans l'espace reel (cosddh, ddelta) pour ; une calibration directe, ; . Cygne dans l'espace interferometrique pour une ; simulation dans une calibration auto. pix_obj, $ ; pixel objet (radians), definie dans SIMUL_VISIB_1 visib_i, $ ; TF (centree dans espace TF) de l'objet simul'e en po- ; sition reelle calculee sur la grille carree conju- ; guee de celle sur laquelle l'objet est defini. ; Dans le cas de l'appel par VISIB_CYG pour une ; calibration directe, le modele spatial est defini ; dans l'espace reel (cosddh, ddelta). ; visib_i a la meme dimension (512*512) que objet_i ; puisqu'elle en est la TF par fft. pix_uv ; valeur du pixel (radians^-1) sur visib. stop end ; fin de T_SIMUL_VISIB_1 ;______________________________________________________________________________ PRO SIMUL_VISIB_1, $ ; simulation de source pour CALIBRATION. ; entree i_calib, $ ; 1 pour Cygne calib directe, 2 pour pseudo Cygne. i_champ_uv , $ ; usuel 1, 2 pour tracer la TF -> 2eme maximum ; (appel par VIS_MOD_CYG (dans fichier du meme nom)). ; sortie objet_i, $ ; objet simul'e en position reelle et defini sur usuel- ; ment 512*512 pts sur un champ 22' (2048*2048 pts ; sur un champ 88' d'arc jusqu'au 8 nov 04) pour ; etre bien decrit : ; . Cygne dans l'espace reel (cosddh, ddelta) pour ; une calibration directe, ; . Cygne dans l'espace interferometrique pour une ; calibration auto. pix_obj, $ ; valeur du pixel (radians) de l'objet en sortie, uti- ; lisee seulement dans VISIB_CYG (dans CALIB_SUB.PRO) visib_i, $ ; Complexar(2048, 2048) : visibilite sur la grille ; conjuguee de la grille com- ; pletee des harmoniques du RH. ; Pour une calibration directe, SIMUL_VISIB_1 est appe- ; le par VISIB_CYG (fichier CALIB_SUB.PRO). Le mode- ; le spatial est defini dans l'espace reel (cosddh, ; ddelta). ; visib_i est la TF (centree dans espace TF) calculee ; sur la grille carree conjuguee de celle sur laquel ; le l'objet est defini dans l'espace direct, et ; dont le pixel a ete choisi de 2e-4 rad. ; Pour une calibration auto SIMUL_VISIB_1 est appele ; (par SIMUL_HARM_N, lui meme appel'e par ; RH_DPATCHFITS_NRH. ; Visib est defini dans l'espace interferometrique ; Pour une calibration de type "auto". pix_uv ; valeur du pixel (radians^-1) sur visib_i. ; ; But : obtenir la visibilite du Cygne pour une calibration : ; A) calibration directe : ; . SIMUL_VISIB_1 est toujours appele par VISIB_CYG pour obtenir un mo- ; dele de visibilite a comparer aux donnees. ; . si les donnees sont simules SIMUL_VISIB_1 est appel'e au prealable ; par SIMUL_HARM_N, lui meme appel'e par RH_DPATCHFITS_NRH. ; B) calibration de type auto : SIMUL_VISIB_1 est appele (par SIMUL_HARM_N ; lui meme appel'e par RH_DPATCHFITS_NRH) que si les donnees sont ; simulees. On peut dans ce cas ajouter une base large au Cygne. ; Usage courant (depuis la mise en service de la nouvelle calibration en fev ; 2003) : etre appele par VISIB_CYG, lui meme appele par CALIBRATION ; pour faire une calibration directe ; ; Dans tous les cas la visibilite simulee est celle de l'objet source : ; - pour une simulation directe les valeurs de la visibilite a comparer aux ; donnees (reelles ou simulees) sont ensuite calculees par VISIB_CYG, ; - pour une simulation de donnees (en vue d'une calibration directe ou simu- ; lee) le vecteur harm_N simul'e est produit par SIMUL_HARM_N, qui lui ; meme appelle VISIB_CYG s'il s'agit d'une simulation de donnees pour une ; calibration directe. ; Etapes de SIMUL_VISIB_1 en usage courant (redige le 20 nov 2003) ; Choix : ; dimension d'image : np = 128 pts, puis interpolation a 512. ; pixel 5 e-5 rad = 10.3" => champ = 22.0' ; definition source dans espace reel ; ; calcul modele : ; calcul d'une composante interpolee facteur f_i pour immersion dans ; tableau (f_i*np, f_i*np) ; calcul des 2 composantes dans le champ (f_i*np, f_i*np), avec barycen- ; tre au centre du champ de np points interpole par facteur f_i=4. ; translation du barycentre en position reelle (en fait translation nulle ; pour une calibration directe). ; calcul de la visibilite avec origine au centre du champ (a ce stade le ; champ TF est de np*f_i = 512 points, et la visibilite s'etend bien ; au dela de la bande du RH) ; En sortie visib est donc la visibilite sur une grille a mailles carrees, ; conjuguee de la grille de l'espace direct de pixel 5e-5 rad. ; Relais par VISIB_CYG (fichier CALIB_SUB) : "visib" est transform'e en ; "visib_RH", visibilite (centree en milieu de tableau) sur la grille ; completee des harmoniques du RH. ; "visib_RH" inclut la decroissance en freq^-0.9 alors que "visib_i" ne ; dans SIMUL_VISIB_1 l'inclut pas. ; ; Modifications: ; 02 mar 26: introduction de i_calib dans la liste d'entree ; i_calib=1 pour une calibration directe. ; on determine un modele realiste des harmoniques du Cygne ; dans le plan des frequences spatiales conjuguees de cosd*dh ; et ddelta. On choisit le pixel (en radians) pour decrire ; l'objet CYGNE avec plusieurs pixels. ; Remarques : ; . visib_i est la TF du Cygne sur une grille carree. Le calcul ; de cette TF sur une grille contenant les harmoniques du ; RH est faite dans VISIB_CYG. ; . ca conviendrait aussi a une calibration de type auto. ; i_calib=2, tolerable mais non necessaire pour une calibration ; de type auto. ; La source est decrite en pixels interferometriques, sans re- ; ference a la position (H delta) de la source ni a la fre- ; quence d'observation. ; Cette option est appelee a disparaitre mais a l'avantage de ; preserver la 1ere utilisation en simulation. ; 02 sep 27 : on calcule la TF sur npi (2048) points et non plus sur np (128) ; points pour eviter des artefacts dus au fait que le modele de ; source est alors defini sur trop peu de points. ; 03 sep 3 : reglage des parametres definnssant le Cygne, par comparaison avec ; l'observation du 26 aug 03. Le detail de la comparaison est ; decrit dans Cygne.txt. ; 04 nov 8 : . la valeur du pixel pour le calcul de l'objet etait rest'e egal a ; 0.0002 rad (41" d'arc). On le reduit a 5 10^-5 rad (10.3"). ; La valeur du coeff d'interpolation f_i (pour bien decrire ; l'objet "Cygne", avec une distance entre les composantes ; arrondie a + ou - 1" d'arc), peut ainsi etre ramenee de 16 ; a 4, ce qui fait une TF sur 512 points, plus rapide que sur ; 2048. ; . correction de la faute trouvee : la reduction de la TF a np*np ; points (pour avoir la grille completee des harmoniques du RH) ; etait incorrectement faite en ne gardant qu'un pave central ; dans la TF 2048*2048. Il faut en fait calculer les coordonnees ; des points ou il faut "piquer" les freq spatiales du RH. ; C'etait en fait ce qui etait fait dans VISIB_CYG, seul le com- ; mentaire (dans SIMUL_VISIB_1) que cette reduction donnait la ; grille du RH etait faux. Ce n'etait pas trop grave car il se ; trouvait qu'on avait pris un pixel de 0.0002 rad = 41" d'arc, ; et donc un champ de ; 41" * 128 = 88' d'arc = 0.0256 rad, ; ce qui serait le champ a 224 MHz. L'extension le la TF reduite ; etait donc suffisante pour que VISIB_CYG en extrait la grille ; des harmoniques du RH. ; Ca serait devenu tres faux avec le nouveau pixel de 10.3". ; => comme SIMUL_VISIB_1 ne peut pas calculer la grille des ; harmoniques du RH car il n'a ni l'heure ni la declinaison, on ; reporte de travail d'extraction en aval, dans VISIB_CYG. ; VISIB_CYG est appele par : ; . CALIBRATION pour la calibration, ; . SIMUL_HARM_N pour la simulation. ; On laisse un peu tomber le cas d'une calibration "auto" (avec ; icalib=2). ; . modif de la liste de sortie (np en moins, pix_uv en plus) ; if (i_calib eq 2) then print, $ ; 'SIMUL_VISIB_1 : simulation de source pour une autocalibration' ; Distinction entre les appels par VISIB_CYG et SIMUL_HARM_N ; on definit la valeur du pixel pix_xy, dans l'espace direct. ; rappel : . i_calib=1 les dimensions sont definies en secondes d'arc, ; . i_calib=2 ------------------------------- pixels. ; Pour i_calib=1 on convient de prendre ; pix_xy = 5 10^-5 rad = 10.31", ; ce qui donne un champ (sur 128 points) de ; champ = 0.0064 rad = 22.0' d'arc ; suffisant pour decrire le Cygne. ; En fait pour bien decrire le modele de Cygne, ce modele est interpole ; d'un facteur f_i = 4 (512 * 512 points sur le champ de 22'd'arc) ; ; La TF (donc aussi 512 * 512) est calculee sur un grille de pas : ; pix_uv = 1/champ = 1/(np * pix_xy) = 156 rad^-1 ; La freq spatiale max est : ; np/2 * pix_uv = np/2 / (np*pix_xy) = 1 / (2*pix_xy) ; soit 10 000 rad^-1 . ; Choix de la dimension des tableaux et du pixel objet np = 128 ; nbre de points pour decrire l'objet> Sera interpole ; par le factreur f_i. Je conserve les 2 etapes pour ; des raisons historiques. if (i_calib eq 1) then begin pix_xy = 0.00005 ; pixel (rad) soit 10.3" d'arc. pix_xy_s = pix_xy * 3600 * 180/!pi ; -------- secondes d'arc. endif if (i_calib eq 2) then begin pix_xy = 1 ; bidon. endif champ_obj = np * pix_xy ; radians. pix_uv = 1 / champ_obj ; radians^-1 ; Definition du modele de source et de sa position ; Composante compacte ; les dimensions doivent etre comprises : ; - en pixels (donc inchangees) si i_calib=2 (appel par SIMUL_HARM_N) ; - en minutes d'arc si i_calib=1 (appel par VISIB_CYG). Dans ce cas ; la valeur du pixel est 5 10^-5 radian = 0.172' = 10.3" d'arc. xg = np/2 ; abscisse centre gravite (entier, 0 a np-1). yg = np/2 ; ordonnee -------------------------------- if (i_calib eq 1) then begin ; simuler Cygne dans espace reel. ; Valeurs adoptees le 30 jan 87 (calibration, calculs litteraux de 85-87) : ; dd = 81.27 ; distance entre composantes (sec d'arc). ; alpha = 23.75 ; inclinaison ligne des sources (deg). ; rap = 0.812 ; rapport intensites droite/gauche. Usuel=0.8 . ; larg_c = 8.51 ; 1/2 larg des sources a mi-hauteur (" d'arc). ; Valeurs adoptees le 3 sep 03, a partir du Cygne du 26 aug 03 avec compa- ; raison avec passage du VIS_MOD_CYGNE ; (fichier du meme nom). dd = 82 ; distance entre composantes (sec d'arc). alpha = 19.5 ; inclinaison ligne des sources (deg). rap = 0.80 ; rapport intensites droite/gauche. Usuel=0.8 . larg_c = 16.0 ; 1/2 larg des sources a mi-hauteur (" d'arc). ; Conversion en pixels. dd = dd / pix_xy_s larg_c = larg_c / pix_xy_s endif if (i_calib eq 2) then begin; simuler Cygne dans espace interferometrique. dd = 2.01 ; distance entre composantes (pixels, sauf 0). alpha = 24 ; inclunaison ligne des sources (deg). Usuel 24 rap = 0.7 ; rapport intensites droite/gauche. usuel=0.8 . larg_c = 1.0 ; 1/2 larg des sources a mi-hauteur (pixels). ; parametres pour essais: dd = 1.50 rap = 0.5 larg_c = 1.3 ; 0.01 2.0 endif flux_c = 1.0 ; flux total (SFU). Pour le Cygne on tient ; compte de freq^-0.9 dans VISIB_CYGNE. alpha = !pi / 180 * alpha ; Composante large xb = 64 ; abscisse centre gravite base (entier, 0 a np-1). yb = 64 ; ordonnee ------------------------------------- larg_ew = 40 ; 1/2 largeur EW hors tout en pixels. larg_ns = 35 ; ----------- NS ------------------- ae = 3 ; exposant ( fonction (1 - x^ae)^ae ). Usuel 3 ou 4. if (i_calib eq 1) then begin flux_b = 0.0 ; pas de base pour le vrai Cygne. endif if (i_calib eq 2) then begin flux_b = 0.01 ; flux total (SFU). Ne pas depasser 30 fois flux_c ; (pb avec calcul largeur mi-puissance de la source ; compacte au dela) endif ; Remarque : meme a 169 l'image est un peu aliasee en EW au meridien. ; il faut donc prendre une 1/2 largeur totale > 32 pixels. Avec le ; modele (1-x^4)^4 on a : ; 1/2 aliasing 0.02 0.05 0.07 0.10 0.15 0.20 ; 1/2 larg (pixels) 36 37 39 40 41 42 ; A 327 et 435 MHz l'image est tres aliasee en EW (1/2 larg-ew ~ 50 ; pixels). ; Largeur NS : le champ angulaire NS du RH au meridien est proportion- ; nel a : ; cos(teta) / Dns ; ou teta est la distance zenithale. ; delta=23 delta=0 delta=-23 ; cos(teta)=0.91 cos(teta)=0.68 cos(teta)=0.34 ; cos(teta) est donc toujours < Dew / Dns et au meridien le champ NS ; est toujours au champ EW avec E0. Le diametre solaire NS etant deja ; inferieur a som diametre EW on doit donc prendre larg_ns nettement ; plus petit que larg_ew. Mais des des qu'on s'ecarte du meridien le ; champ NS augment. En moyenne on prend larg_ns un peu plus petit. ; Calcul des positions (en pixels par rapport au centre de gravite dans l'objet ; np*np) et des amplitudes relatives des composantes de la source ; compacte. d1 = - (rap *dd) / (rap + 1) d2 = dd / (rap + 1) x1 = d1 * cos(alpha) ; en pixels de l'image sur np pts. y1 = d1 * sin(alpha) x2 = d2 * cos(alpha) y2 = d2 * sin(alpha) a1 = 1 / (1 + rap) ; amplitude relative 1ere composante. a2 = rap / (1 + rap) ; ------------------ 2eme ----------. ; Calcul de la distribution de brillance de la source compacte ; Se fait en plusieurs etapes : ; - calcul d'une composante situee au centre d'un champ reduit avec une ; interpolation (facteur f_i) pour pouvoir faire des translations ; precises et immersion dans un grand tableau (f_i*np, f_i*np). ; - calcul des 2 composantes de la source compacte avec la meme inter- ; polation ; - reduction son centre ; de gravite au milieu d'un champ carre comprenant np1 pixels de part ; et d'autre de son centre, c'est a dire s'etendant sur une grille de ; (2*np1 + 1) * (2*np1 + 1) pixels. On prends le centre de ce champ ; reduit en (np/2, np/2) et on choisit np1 assez petit pour eviter les ; underflows dans le calcul des exponentielles. ; - pair (precaution du 9 oct 00 pour eviter des decalages intempes- ; tifs lors de l'immersion) ; On decrit la source compacte par 128*128 pts sur ce champ reduit (in- ; terpolation par un facteur 128/np1). ; Calcul de la composante elementaire de la source compacte dans un champ ; reduit interpol'e np1 = nint(0.0*dd + 6.0*larg_c) ; np1 1/2 largeur (en canaux initiaux) ; du champ reduit. 6 *larg_c pour ; ne pas couper les pieds de la ; gaussiennes sans aller trop loin ; (underflows). 0.0*dd historique. if (np1 lt 2 ) then np1 = 2 ; champ reduit entier, force a etre ; non nul. if (np1 gt np/4) then np1 = np/4 f_i = 4 ; facteur d'interpolation. np3 = np1 * f_i ; 2*np3 + 1 nbre de points sur le ; champ reduit interpole. pix_obj = pix_xy / f_i ; pixel sur l'objet en sortie. d_c = shift(dist(2*np3 + 1), np3, np3) / f_i ; d_c distance au centre dans le champ reduit interpole, expri- ; mee en pixels initiaux. Le pas de variation sur les axes, ; en pixels initiaux,est 1/f_i. sig = larg_c / sqrt(2 * alog(2)) ; ecart type de la gaussienne ; exp(-x^2 / (2 * sig^2)) arg = - d_c^2 / (2 * sig^2) ; argument gaussienne. dom_u = where (arg lt -20, ic) if (ic gt 0) then arg(dom_u) = -20 ; pour eviter les underflows. amp = exp(arg) ; composante elementaire, decrite par ; (2*np1 + 1) * (2*np1 + 1) pts sur ; le champ reduit. ; Controle du bon positionnement de amp icon = 0 if (icon eq 1) then begin print, 'SIMUL_VISIB_1 : controle du bon positionnement de "amp"' min_arg = min(arg) help, pix_xy_s, dd, larg_c, np1, np3, sig, arg, min_arg, ic absc = (- np3 + indgen(2 * np3 + 1)) / float(f_i) ; en canaux initiaux. window, 0 shade_surf, arg xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Argument de l''exponentielle' window, 1 plot, absc, amp(*, np3) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Coupe EW d une composante de source compacte' window, 2 plot, absc, amp(np3, *) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Coupe NS d une composante de source compacte' stop ; Conclusion (6 oct 2000) : le max est bien en 0, cad au milieu. endif ; Immersion de amp dans un champ complet interpol'e ; la largeur de ce champ interpol'e est toujours pix_xy * np (en fait ; presque pix_xy * (np+1) du fait de l'interpolation) npi = np * f_i ; Attention ici npi n'est pas le nbre de pts sur l'image interferome- ; trique (128). amp_i = fltarr(npi, npi) amp_i(npi/2 - np3, npi/2 - np3) = amp ; rappel : np3 1/2 larg du champ reduit (sur lequel on calcule une ; composante) interpole par le facteur f_i ; Controle de amp_i icon = 0 if (icon eq 1) then begin print, ' de amp dans champ complet interpol''e' toto = amp_i(npi/2-np3:npi/2+np3, npi/2-np3:npi/2+np3) absc = (- np3 + indgen(2 * np3 + 1)) / float(f_i) ; en canaux initiaux. window, 0 & wset, 0 plot, absc, toto(*, np3), xstyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Coupe EW source elementaire immergee (partie centrale) window, 1 & wset, 1 plot, absc, toto(np3, *), xstyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Coupe NS source elementaire immergee (partie centrale) stop endif ; Calcul source compacte au centre du champ complet interpole x1_i = nint(x1 * f_i) ; absc 1ere composante, pixels image interpolee y1_i = nint(y1 * f_i) x2_i = nint(x2 * f_i) y2_i = nint(y2 * f_i) obj_c_c_i = a1 * shift(amp_i, x1_i, y1_i) + a2 * shift(amp_i, x2_i, y2_i) ; comme objet compact au centre interpole. amp_i = 0 ; destruction de amp_i (economie de place). ; Normalisation de la composante compacte (30 sept 02) ; La normalisation a pour but que les valeurs de la TF soient correctes ; (composante continue et bas harmoniques egaux au flux en sfu). ; Or la TF sera calculee sur npi = np * f_i points (np=128 et f_i=16) ; de facon qur l'objet soit decrit par un nombre de points raisonna- ; ble, puis on ne gardera que le carre central 128*128 du champ TF. ; La somme des points de l'objet defini sur npi points doit donc etre ; flux_c (et non flux_c * f_i^2), puisque en IDL : ; fft inverse a l'origine = Somme des points ; Rappel sur fft inverse : drapeau +1, signe + dans l'exponentielle. ; Si a est le facteur de normalisation : ; total(a * obj_i) = flux_c ; => a = flux_c / total(obj_i) obj_c_c_i = flux_c / total(obj_c_c_i) * obj_c_c_i ; Controle de obj_c_c_i icon = 0 if (icon eq 1) then begin print, 'SIMUL_VISIB_1 : controle de la composante compacte interpolee' print, ' avec barycentre a l''origine' ch_format = "(17x, 'attention : profils EW et NS alors que " + $ "alpha =', f5.1, ' deg.')" print, format=ch_format, 180/!pi * alpha print, ' (obj_c_c_i), traces de profils EW et NS sur +-', np1, $ ' canaux initiaux' absc = (- np3 + indgen(2 * np3 + 1)) / float(f_i) ; en canaux initiaux. ; rappel 2*np3 + 1 nbre de points sur le champ reduit de calcul de ; la source compacte interpolee. window, 0 nc = npi/2 ; pt central du champ total interpole. plot, absc, obj_c_c_i(nc-np3:nc+np3, nc), xstyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Partie centrale coupe EW source compacte complete' window, 1 plot, absc, obj_c_c_i(nc, nc-np3:nc+np3), xstyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Partie centrale coupe NS source compacte complete' print, ' total(obj_c_c_i) =', total(obj_c_c_i) stop ; Conclusion (6 oct 2000) : le max est bien en np*f_i, cad au milieu. ; Verif du 27 sept 02 avec calcul du barycentre : ; xbar = 0.006 ybar = 0.024 ; (ce calcul prend 30 sec) endif ; Calcul source compacte en sa position reelle dans le champ interpole im_c_i = shift(obj_c_c_i, f_i*xg - npi/2, f_i*yg - npi/2) obj_c_c_i = 0 ; destruction de obj_c_c_i (economie de place). ; Rappels : ; . npi = np*f_i nbre de pts sur objet interpol'e. ; . xg et yg sont les coord du barycentre de la source compacte. Ils ; sont comptes a partir du coin inferieur gauche et non a partir ; du centre du champ. ; Quand l'objet etait defini sur np (np=128) points il fallait ; imperativement xg et yg entiers. Ce n'est plus necessaire. ; Saut du calcul de la composante large dans le cas du vrai Cygne (i_calib=1) if (i_calib eq 1) then begin base_i = 0 ; base bidon a_tf_b = 0 ; abs(tf(base)) bidon. goto, a1 endif ; Calcul de la distribution de brillance de la source large. ; On utilise une repartition du type : ; amp = (1 - r^p)^p pour r<1 (1) ; amp = 0 pour r>1 ; ou r est un rayon vecteur normalise (max = 1) et p un exposant. ; Dans une distribution anisotrope elliptique les lignes de niveau ; sont des ellipses homothetiques d'equation ; x^2/a^2 + y^2/b^2 = 1 (2) ; On pose a = m * A et b = m * B ou A et B sont les 1/2 axes de ; l'ellipse la plus grande (de niveau nul). On a donc ; x^2/A^2 + y^2/B^2 = m^2 (3) ; Le 1er membre de (3) est donc constant sur une ellipse donnee et on ; peut le prendre comme "rayon vecteur" r dans l'eq (1) ; Definition du champ reduit carre dans lequel la souce large est centree ; La source etant large par hypothese, on la definit directement sur ; np2*np2 pts, sans interpolation. On la calcule sur un champ reduit, ; centre au centre su soleil calme, et dont la largeur est limitee ; par le bord du champ le plus proche. np2 = 2 * min([np-xb, xb, np-yb, yb]) ; Calcul des coordonnees centrees et reduites, x/A et y/B et du "rayon vec- ; teur r limit'e a [0, 1]. x_r = (findgen(np2) - np2/2) / larg_ew ; x/A notation precedente. y_r = (findgen(np2) - np2/2) / larg_ns ; y/B ------------------- mx = x_r # replicate(1, np2) my = replicate(1, np2) # y_r re = sqrt(mx^2 + my^2) ; Limitations a de x_r, t_r et re dom = where(re gt 1, ic) if (ic gt 0) then re(dom) = 1 amp = (1 - re^ae)^ae ; Immersion dans un tableau np*np. ima_b = fltarr(np, np) ; b comme base. ima_b(np/2-np2/2, np/2-np2/2) = flux_b / total(amp) * amp ; Controle de source large simulee au centre du champ total (np/2, np/2). icon = 0 if (icon eq 1) then begin window, 0 & wset, 0 im_visu = congrid(ima_b, 512, 512, cubic=-0.5) shade_surf, im_visu xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'source large simulee au centre du champ total' window, 1 & wset, 1 tvscl, im_visu stop endif ; Objet source large a sa position finale im_b = shift(ima_b, xb - np/2, yb - np/2) ; b comme base. ; Controle de source large simulee a sa position finale. icon = 0 if (icon eq 1) then begin help, im_b ; Traces 2D window, 0 & wset, 0 shade_surf, im_b xyouts, 0.5, 0.07, /normal, alignment=0.5, $ 'Base large sur tout le champ' window, 1 & wset, 1 tvscl, congrid(im_b, 512, 512, cubic=-0.5) xyouts, 0.5, 0.07, /normal, alignment=0.5, $ 'Base large sur tout le champ' ; Traces 1D absc = - np/2 + indgen(np+1) window, 2 & wset, 2 plot, absc, im_b(* , np/2), xstyle=1 oplot, absc, im_b(np/2, * ), linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Source large position finale : coupe EW (___) et NS(...)' stop endif ; Calcul TF de la base (utilisee plus bas pour limiter domaine de trace) tf_b = shift(fft(shift(im_b, -np/2, -np/2), 1), np/2, np/2) ; tf_b TF quasi-reelle et centree de la base. a_tf_b = abs(tf_b) ; Trace de la TF de la base (pour voir s'il y a une remontee lointaine) icon = 0 if (icon eq 1) then begin print, 'SIMUL_VISIB_1 : trace TF de la base' nz = 16 ch_nz = string(nz) ; zoom atfb = a_tf_b(np/2-nz:np/2+nz, np/2-nz:np/2+nz) atfb_512 = congrid(atfb, 512, 512, cubic=-0.5) window, 0 & wset, 0 shade_surf, atfb_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(TF base) harm ' +ch_nz window, 1 & wset, 1 tvscl, atfb_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(TF base) harm ' +ch_nz window, 2 & wset, 2 absc = -nz + indgen(2*nz + 1) plot, absc, atfb(* , nz), xstyle=1, yrange=[0, 2] oplot, absc, atfb(nz, * ), linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(TF base) : coupes EW (___), et NS (...)' stop ; Conclusion (13 oct 2000) : la visibilite ne remonte pas au dela de 5. endif ; Interpolation et normalisation de la source large sur npi points (usuel 2048) ; en vue d'additionner a la source compacte, definie aussi sur npi pts. ; Ceci pour calculer la TF sur npi pts, de facon a eviter les erreurs ; de discretisation qui conduisaient a un barycentre decale, a une ; pente de phase et une image du Cygne, calcul'e sur la grille (comple- ; tee) des harmoniques du RH, qui est mauvaise a basse frequence ; (essais du 26 sept 02), avec une distorsion et un rapport des compo- ; santes visiblement biais'e. base_i = congrid(im_b, npi, npi, cubic=-0.5) ; Rappel : la base large est aussi a sa popsition finale, definie par ; xb et yb par rapport au coin gauche du champ sur np (128) points). if (total(base_i) gt 0) then base_i = flux_b / total(base_i) * base_i ; En effet lors de la precedente normalisation sur np points im_b peut ; etre nul si flux_b=0. a1: J_M = 0 ; fin du saut du calcul de la composante large si i_calib=1. ; Objet final (base large + composante compacte) : ; - en sa position reelle, ; - avec la source compacte au coin inf gauche, en vue d'obtenir une ; TF quasi-reelle pour une visualisation commode. En fait si la ; source compacte n'est pas au centre du disque, les bas hqrmoni- ; ques, domines par le soleil calme, ne sont plus quasi-reels. ; Modif du 27 sept 02 : on calcule la TF sur npi points et non plus ; sur np points pour eviter des artefacts dus au fait que le ; modele de source est alors defini sur trop peu de points (voir ; plus haut) objet_i = base_i + im_c_i ; respectivement base large et source compacte. ; rappel : si i_calib=1 base_i = 0 im_c_i = 0 ; destruction de im_c_i (economie de place). base_i = 0 ; -------------- base_i ------------------- ; Controle objet final icon = 0 if (icon eq 1) then begin print, 'SIMUL_VISIB_1, controle objet final interpol''e, pos. reelle' print, ' rappel : facteur d''interpolation =', f_i nz = 10 ; zoom sur -nz a +nz pixels initiaux. objet_red = objet_i(npi/2 - nz*f_i : npi/2 + nz*f_i, $ npi/2 - nz*f_i : npi/2 + nz*f_i) window, 0 im_visu = congrid(objet_red, 512, 512, cubic=-0.5) shade_surf, im_visu ; , az=-20 ; rotation z champ_red = 2 * nz * pix_xy_s / 60 ; champ en min d'arc. ch_format = "(f5.2, ' min d arc')" ch_champ = string (format=ch_format, champ_red) xyouts, /normal, 0.5, 0.97, alignment=0.5, $ 'modele final , zoom sur un champ de : ' + ch_champ xyouts, /normal, 0.5, 0.92, alignment=0.5, $ '(position reelle)' window, 1, xsize=512, ysize=512 tvscl, im_visu xyouts, /normal, 0.5, 0.97, alignment=0.5, $ 'source finale, zoom sur un champ de : ' + ch_champ xyouts, /normal, 0.5, 0.92, alignment=0.5, $ '(position reelle)' ; window, 2 & wset, 2 ; im_visu = congrid(objet_coin, 512, 512, cubic=-0.5) ; shade_surf, im_visu ; xyouts, /normal, 0.5, 0.95, alignment=0.5, $ ; 'source totale avec source compacte au coin inf gauche' ; window, 3 & wset, 3 ; tvscl, im_visu ; xyouts, /normal, 0.5, 0.95, alignment=0.5, $ ; 'source totale avec source compacte au coin inf gauche' objet_red = 0 ; destruction de objet_red (economie de place). help, objet_i print, ' total(objet final interpol''e) =', total(objet_i) stop endif ; Calcul de TF d'objet source en position finale et de sa visibilite centree. ; La TF se fait ici sur l'objet bien decrit : ; . pixel initial = 0.00005 rad = 10.3" a partir du 8 nov 04 (au lieu ; de 41" d'arc avant le 8 nov 04). ; . on a interpole d'un facteur 16, ce qui reduit le pixel interpole ; a 0.63 " d'arc. ; ; Rappel sur les fft en IDL : ; - les origines dans les espaces direct et de Fourier ne sont pas ; au centre des tableaux mais a leur debut. ; - la TF directe (drapeau -1) comporte un facteur 1/n et un ; signe - dans l'exponentielle complexe. ; - la TF inverse (drapeau 1) ne comporte pas de facteur 1/n et un ; signe + dans l'exponentielle complexe. ; - Proprietes : ; . TF inverse(0) = somme des points de la fonction de depart. ; . total(TF directe) = fonction de depart en zero ; . TF directe(TF invserse) = TF invserse (TF directe) = identite ; . TF directe(TF directe) = 1/n * identite ; . TF invserse(TF invserse) = n * identite ; => la TF inverse realise l'operation faite par le recepteur RH, avec ; une composante continue egale au flux (sfu) et des harmoniques ; egaux au flux de la source si elle est ponctuelle. ; la TF directe (drapeau -1) calcule l'image de la source qui donne ; de tels harmoniques. C'est ce qui est fait dans MALC_IN_2D : ; image = float (shift(fft(harm, -1), np/2, np/2)) ; Rappel : depuis le 27 sept 02 on fait la TF sur un tableau de npi ; points (npi=2048) et non plus sur un tableau de np points ; (np=128) ch_npi = strcompress(string(format="(i4)", npi), /remove_all) ; Rappel : /remove_all equivaut a /remove=all . i_imp_temps_calcul = 0 if (i_imp_temps_calcul eq 1) then $ print, 'SIMUL_VISIB_1 : debut TF sur ' + ch_npi + ' * ' + ch_npi + $ ' points (objet interpole -> TF)' t0 = systime(1) ; temps en sec depuis le 1er jan 70. visib_i = shift( fft( shift(objet_i, -npi/2, -npi/2), 1), npi/2, npi/2) ; role des parentheses fft 1er shift ; environ 8 sec sur mesopz pour calculer une TF 2048*2048. t1 = systime(1) ch_format ="(' fin --------------------------------" + $ "---------- (', f5.2, ' sec)')" if (i_imp_temps_calcul eq 1) then print, format=ch_format, t1 - t0 icon = 0 ; traces shade_surf et tvscl de la visibilite. if (icon eq 1) then begin window, 0 shade_surf, congrid(abs(visib_i), 512, 512, cubic=-0.5) print, ' total(objet_i) =', total(objet_i) print, ' visib_i(npi/2, npi/2) =', float(visib_i(npi/2, npi/2)) window, 1, xsize=512, ysize=512 tvscl, congrid(abs(visib_i), 512, 512, cubic=-0.5) print, ' total(objet_i) =', total(objet_i) print, ' visib_i(npi/2, npi/2) =', float(visib_i(npi/2, npi/2)) help, objet_i, visib_i stop endif ; Ancienne place de la reduction (sans justifications) de la TF a un champ de ; np*np points. le calcul de la TF sur la grille completee des harmo- ; niques du RH est reportee dans : ; . VISIB_CYG pour la calibration ; . SIMUL_HARM_N pour la simulation. ; visib = visib_i ; visib(0, 0) = visib_i(npi/2 - np/2 : npi/2 + np/2 - 1 , $ ; npi/2 - np/2 : npi/2 + np/2 - 1 ) ; supprime le 8 nov 04. ; De meme, on supprime le 8 nov 04 la reduction pour l'usage special a ; VIS_MOD_CYG, ou il ne s'agit que de visualiser la visibilite. ; Reduction a un champ de TF de 2np*2np points pour application speciale a ; appel par VIS_MOD_CYG (dans fichier du meme nom): il faut aller ; assez loin sur la TF pour observer le 1er maximum lateral (on ; cherche a ajuster le parametre larg_c qui gouverne la decrois- ; sance lente) ; if(i_champ_uv ne 1) then begin ; i_champ_uv=2 dans le cas de l'appel ; ; par VIS_MOD_CYG (dans fichier du ; ; meme nom). ; npuv = i_champ_uv * np ; visib = complexarr(npuv, npuv) ; visib(0, 0) = visib_i(npi/2 - npuv/2 : npi/2 + npuv/2 - 1 , $ ; npi/2 - npuv/2 : npi/2 + npuv/2 - 1 ) ; endif ; La suite ne consiste qu'en controles. ; Controle sommaire pour le cas utilisation a partir de VIS_MOD_CYGNE (dans ; fichier du meme nom). icon = 0 if (icon eq 1) then begin window, 0 shade_surf, abs(visib_i) help, visib stop endif tf_obj_coin_c = visib_i ; TF centree. Rustine pour faire tourner les controles ci-dessous. ; En fait, ca n'est pratique que si tf_obj_coin_c est la TF qusai- ; reelle de 'objet centre, alors que visib est la TF de l'objet en ; position reelle. => prendre xg = yg = xb = yb = 64 pour faire tourner ; ces controles. ; Definition du domaine d_ph accepte pour le trace de la phase de la source ; totale (qui sera mise a zero ailleurs) de facon a eviter les sauts ; de 360 degres dus a la visibilite de la base ; On verifie au moins une des 2 conditions : ; - visibilite > fraction du flux de la base ; - ---------- > quelques fois la visibilite de la base. ; la 1ere condition intervient pres de l'origine quand la base est ; forte ou seule, la 2eme condition intervient pour les harmoniques ; de la source compacte surpassent ceux de la base. a_vis = abs(visib_i) af1= 0.6 ; limite d'amplitude pour eviter les sauts de phase, ; qui apparaissent pour 0.05 avec ae=3 et flux_c=0. af2 = 2 ; pour que la source compacte surpasse la base. d_ph = where((a_vis gt af1*flux_b) or (a_vis gt af2*a_tf_b), ic_d_ph) d_ph = [d_ph, np/2*np + np/2] ; rajouter (np/2, np/2) ; Controle de d_ph icon = 0 if (icon eq 1) then begin toto = replicate(0, np, np) if (ic_d_ph gt 0) then toto(d_ph) = 1 window, 0 & wset, 0 shade_surf, toto xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Domaine pour tracer phase de source totale' stop endif ; Controle de TF (recentree) avec la source compacte au coin inferieur gauche. ; Traces a 2 dimensions icon = 0 if (icon eq 1) then begin print, 'SIMUL_VISIB_1 : controle de la visibilite de la source' nz =60 ; nbre d'harmoniques a zoomer. nd = np/2 - nz & nf = np/2 + nz ; alleger notations. ch_nz = string(nz) tf_z = tf_obj_coin_c ; On ne zoome pas a ce stade pour pouvoir ensuite mettre la phase ; a zero sur le domaine d_ph. amp = abs(tf_z) phas = 180/!pi * atan(imaginary(tf_z), float(tf_z + 0.0001)) ; 0.0001 pour atan. if (ic_d_ph gt 0) then begin amul = replicate(0, np, np) amul(d_ph) = 1 phas = phas * amul endif amp_512 = congrid(amp (nd:nf, nd:nf), 512, 512, cubic=-0.5) phas_512 = congrid(phas(nd:nf, nd:nf), 512, 512, cubic=-0.5) window, 0 & wset, 0 shade_surf, amp_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'abs (TF source) ' + ch_nz + ' harm' window, 1 & wset, 1 tvscl, amp_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'abs (TF source)' + ch_nz + ' harm' window, 2 & wset, 2 shade_surf, phas_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'phase (TF source quasi reelle) deg ' + ch_nz + ' harm' window, 3 & wset, 3 tvscl, phas_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'phase (TF source quasi reelle) deg ' + ch_nz + ' harm' stop endif ; Coupes EW et NS pures icon = 0 if (icon eq 1) then begin window, 0 & wset, 0 nh3 = 12 ; traces harm (-nh3, nh3) absc = indgen(2*nh3+1) - nh3 ; domaine des abscisses. xmin = min(absc) & xmax = max(absc) imin = xmin + np/2 & imax = xmax + np/2 ymax = max(abs(tf_obj_coin_c)) plot, xrange=[xmin, xmax], xstyle=1, $ yrange=[-0.4*ymax, 1.1*ymax], ystyle=1, absc, $ abs (tf_obj_coin_c(imin:imax, np/2)) oplot, absc, float (tf_obj_coin_c(imin:imax, np/2)), psym=2 oplot, absc, imaginary(tf_obj_coin_c(imin:imax, np/2)), psym=3 ; psym=2 *, 3 ., 4 losange xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'TF recentree, harm EW purs : abs ___, reel ***, imag ...' window, 1 & wset, 1 nh3 = 36 ; traces harm (-nh3, nh3) plot, xrange=[xmin, xmax], xstyle=1,$ yrange=[-0.2*ymax, 1.1*ymax], ystyle=1, absc, $ abs (tf_obj_coin_c(np/2, imin:imax)) oplot, absc, float (tf_obj_coin_c(np/2, imin:imax)), psym=2 oplot, absc, imaginary(tf_obj_coin_c(np/2, imin:imax)), psym=3 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'TF recentree, harm NS purs : abs ___, reel ***, imag ...' stop ; tf_obj_coin_c est quasi-reelle (a l'asymetrie de la source pres). endif ; Controle de TF (recentree) de la source en position reelle. Traces 1D. icon = 0 if (icon eq 1) then begin nh = 12 absc = indgen(2 * nh + 1) - nh a_ew_vis = abs (visib(np/2-nh : np/2+nh , np/2)) r_ew_vis = float (visib(np/2-nh : np/2+nh , np/2)) i_ew_vis = imaginary(visib(np/2-nh : np/2+nh , np/2)) a_ns_vis = abs (visib(np/2 , np/2-nh : np/2+nh)) r_ns_vis = float (visib(np/2 , np/2-nh : np/2+nh)) i_ns_vis = imaginary(visib(np/2 , np/2-nh : np/2+nh)) ; Traces coupe EW window, 0 & wset, 0 ymax = max(a_ew_vis) plot, xstyle=1, yrange=[-0.3*ymax, 1.1*ymax], ystyle=1, $ a_ew_vis oplot, r_ew_vis, linestyle=2 oplot, i_ew_vis, linestyle=1 ; oplot, absc, imaginary(tf_obj_coin_c(np/2, np/2-nh:np/2+nh)), psym=3 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'visib simulee, coupe EW, abs __, reel - -, imagin ...' window, 1 & wset, 1 ymax = max(a_ns_vis) plot, xstyle=1, yrange=[-0.3*ymax, 1.1*ymax], ystyle=1, $ a_ns_vis oplot, r_ns_vis, linestyle=2 oplot, i_ns_vis, linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'visib simulee, coupe NS, abs __, reel - -, imagin ...' stop endif ; Controle de TF (recentree) de la source en position reelle. Traces 2D. ; Rappel : a ce stade visib est la TF complete de la source definie ; par le maillage dans l'espace objet. Ce n'est pas encore la TF ; incomplete, definie sur la grille (completee) des harmoniques du ; RH, qui sera calculee par interpolation dans VISIB_CYG. icon = 0 if (icon eq 1) then begin nh = 10 ; pour zoomer autour de l'origine. ch_format = "('(-+ ', i3, ' harm.)')" ch_harm = string(format=ch_format, nh) h1 = np/2-nh & h2 = np/2+nh visib_a = visib(h1:h2, h1:h2) + 0.0001 ; 0.0001 pour atan. a_visib = abs(visib_a) p_visib = 180/!pi * atan(imaginary(visib_a), float(visib_a)) window, 0 shade_surf, a_visib ; $ xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'abs(TF source pos reelle) ' + ch_harm window, 1 shade_surf, az=10, congrid(p_visib, 512, 512, cubic=-0.5) ; az rotation directe autour de l'axe z (defaut -30). xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'phase(TF source pos reelle) ' + ch_harm ; window, 2 ; tvscl, congrid(a_visib, 512, 512, cubic=-0.5) ; xyouts, /normal, 0.5, 0.95, alignment=0.5, $ ; 'abs(TF recentree) source pos reelle' stop ; Rem: visib_a est quasi-reelle pour xg = yg = 64. endif ; Controle de l'effet sur l'objet simul'e de la reduction du champ de la TF ; (de npi a np points) icon = 0 ; 1 si on n'a pas detruit objet_i if (icon eq 1) then begin print, 'SIMUL_HARM_1 : controle de l''effet de limitation de la TF'+$ ' a 128 pts' nz = 5 ; zoom sur -nz a +nz canaux initiaux. image = float(shift(fft(shift(visib, -np/2, -np/2), -1), np/2, np/2)) image_i = congrid(image, npi, npi, cubic=-0.5) ; la TF est assez lisse et on veut voir l'effet de la limitation de ; son domaine (il y a une marche au bord). Un congrid suffit. image_red = image_i(npi/2 - nz*f_i : npi/2 + nz*f_i, $ npi/2 - nz*f_i : npi/2 + nz*f_i) objet_red = objet_i(npi/2 - nz*f_i : npi/2 + nz*f_i, $ npi/2 - nz*f_i : npi/2 + nz*f_i) window, 0 im_visu = congrid(objet_red, 512, 512, cubic=-0.5) shade_surf, im_visu * f_i^2 ; pour comparer avec trace suivant. ch_format = "(i3)" ch_lab = string (format=ch_format, nz) xyouts, /normal, 0.5, 0.97, alignment=0.5, $ 'objet exact , zoom sur - a + ' + ch_lab + ' canaux de base' xyouts, /normal, 0.5, 0.92, alignment=0.5, $ '(position reelle)' window, 1 im_visu = congrid(image_red, 512, 512, cubic=-0.5) shade_surf, im_visu xyouts, /normal, 0.5, 0.97, alignment=0.5, $ 'objet a TF limitee a 128 points; zoom idem' ; window, 2 ; im_visu = congrid(image_red, 512, 512, cubic=-0.5) - $ ; congrid(objet_red, 512, 512, cubic=-0.5) ; shade_surf, im_visu ; xyouts, /normal, 0.5, 0.97, alignment=0.5, $ ; 'objet a TF limitee - objet exact zoom idem' print, ' flux de l''objet :' print, ' exact flux =', $ total(objet_i) print, ' TF limitee a 128 points (sfu) flux =', $ total(image_i)/f_i^2 stop endif end ; fin de SIMUL_VISIB_1 ;______________________________________________________________________________ PRO SIMUL_VISIB_2, $ ; simulation de source pour CALC_IMAGE_HELIO. ; sortie npi, $ ; dimension de l'image (choisir ici le npi de dpnew). soleil, $ ; objet simule en position relle visib ; TF (centree dans espace TF) d'image simulee en ; position reelle. ; But : Simuler la visibilite d'un soleil complexe destinee au test de resti- ; tution de composante continue dans RH_MALC_IM_2D, au test de clean dans ; CALC_IMAGE_HELIO etc. . Le soleil simule consiste en ; . une base de soleil calme, ; . un couloir de filament ; . des sources compactes et moins compactes (toutes engendrees avec des ; sources elementaires identiques, gaussiennes et de largeur larg_c (u- ; suel 1.2 pixels), ; . une arche de CME. ; ; La visibilite simulee est celle de l'objet source, directement defini dans ; l'espace interferometrique. En effet le but est de tester le calcul de ; clean et il est commode de comparer directement l'image nettoyee et l'ob- ; jet initial. ; C'est presque realiste pour le soleil au meridien en ete. ; ; print, 'SIMUL_VISIB_2 : simulation de source pour CALC_IMAGE_HELIO' npi = 128 ; nbre de points (identique au nbre de pts npi de ; l'image interferometrique. ; PARAMETRES REGLABLES DU SOLEIL ; Coefficients multiplicateurs (1 ou zero pour essais) amul_b = 1.0001 ; 1.0 ; soleil calme lisse (NON NUL). amul_fil = 0.5 ; 0.20 ; couloir de filament (evt nul). amul_barre = 1.0 ; 1.0 ; barres allongees (evt nul). amul_lar = 1.0 ; 0.6 ; sources larges (evt nul). amul_ep = 1.0 ; 1.0 ; ------- eparses (evt nul). amul_CME = 0.5 ; 1 pour meme brillance que le SC ; avec intens_CME = 3 ou 4. amul_int = 1 ; ------ intense, 10 pour dominer ; l'image, 100 pour dominer le flux. a_mul_comp = 0.7 ; multiplication d'ensemble des sources ; compactes. Usuel 0.9 ; Autres parametres de modification courante larg_c = 1.2 ; 1/2 largeur a mi-hauteur des sources ; elementaires (pixels). S'applique ; aux composantes de toutes les sour- ; ces (y compris CME) sauf le soleil ; calme. Usuel 1.2, 0.2 pour source ; trs etroite. ; Barres allongees ; Barre plutot horizontale x_bar_h = 52;[52, 52] ; absc centre (pixels) y_bar_h = 68;[68, 52] ; ord -------------- ; Barre plutot verticale x_bar_v = 79;[79, 48] y_bar_v = 59;[59, 59] ; Sources larges xc_lar_1 = 55 & yc_lar_1 = 75 ; coord centre (pixels) xc_lar_2 = 75 & yc_lar_2 = 59 ; -------------------- ; Sources eparses i_simple = 0 ; 0 distribution simple, 1 complete ; Boucle genre CME x_CME = 0.4 ; absc du centre de l'arc pr centre du champ, en unites ; de 1/2 largeur ew du soleil calme (larg_ew). 0.6 y_CME = -0.8 ; ord ------------------------------ (larg_ns) -0.7 r_CME = 15 ; rayon (pixels). dr_CME = 3 ; largeur de l'arc (pixels). puiss = 3 ; degre de la decroissance radiale (usuel 5). intens_CME = 1 ; brillance max du CME en unites du max du SC ; (compte non tenu du raccord). ; Multiplication d'ensemble des sources compactes amul_fil = a_mul_comp * amul_fil amul_barre = a_mul_comp * amul_barre amul_lar = a_mul_comp * amul_lar amul_ep = a_mul_comp * amul_ep amul_int = a_mul_comp * amul_int ; amul_CME = a_mul_comp * amul_CME ; pas vraiment source compacte. ; PARAMETRES FIXES (OU PLUS RAREMENT REGLABLES) DU SOLEIL ; Les flux sont en 0.01 SFU (voir multplication par 0.01 plus bas). ; Definition du soleil calme moyen (suite et fin) xb = 64 ; abscisse centre gravite base (entier, 0 a npi-1). yb = 64 ; ordonnee ------------------------------------- larg_ew = 32 ; 1/2 largeur EW hors tout en pixels. ex : 38, 55 larg_ns = 30 ; ----------- NS ------------------- -- 33, 45 ae = 5 ; exposant ( fonction (1 - x^ae)^ae ). Usuel 3 ou 4. flux_b = 1000. ; flux. Ne pas depasser 30 fois flux_c. ; (pb avec calcul largeur mi-puissance de la source ; compacte au dela) flux_b = amul_b * flux_b ; Remarque : meme a 169 l'image est un peu aliasee en EW au meridien. ; il faut donc prendre une 1/2 largeur totale > 32 pixels. Avec le ; modele (1-x^4)^4 on a : ; 1/2 aliasing 0.02 0.05 0.07 0.10 0.15 0.20 ; 1/2 larg (pixels) 36 37 39 40 41 42 ; A 327 et 435 MHz l'image est tres aliasee en EW (1/2 larg-ew ~ 50 ; pixels). ; Largeur NS : le champ angulaire NS du RH au meridien est proportion- ; nel a : ; cos(teta) / Dns ; ou teta est la distance zenithale. ; delta=23 delta=0 delta=-23 ; cos(teta)=0.91 cos(teta)=0.68 cos(teta)=0.34 ; cos(teta) est donc toujours < Dew / Dns et au meridien le champ NS ; est toujours > au champ EW avec E0. Le diametre solaire NS etant ; inferieur a son diametre EW, on doit donc prendre larg_ns nettement ; plus petit que larg_ew. Mais des des qu'on s'ecarte du meridien le ; champ NS augmente. En moyenne on prend larg_ns un peu plus petit. ; Definition des composantes compactes (abscisses, ordonnees, flux (0.01 SFU), ; largeur) ; Couloir de filament xc_fil = [58, 59, 60, 60, 61, 62, 63, 64, 65, 66, 67, 67, 68, 69, 70] yc_fil = [54, 54, 54, 53, 53, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54] fl_fil = [-3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3] ; Barres allongees (formees d'elements horizontaux et verticaux) ; Element horizontal (5 * 3 pixels) : xc_el_h = [-2, -2, -2, -1, -1, -1, 0, 0, 0, 1, 1, 1, 2, 2, 2] yc_el_h = [-1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1] fl_el = [ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] ; Element vertical (3 * 5 pixels) : xc_el_v = [-1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1] yc_el_v = [-2, -1, 0, 1, 2, -2, -1, 0, 1, 2, -2, -1, 0, 1, 2] fl_el = [ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] xc_1 = [xc_el_h + x_bar_h - 5, xc_el_h + x_bar_h, xc_el_h + x_bar_h + 5] yc_1 = [yc_el_h + y_bar_h - 1, yc_el_h + y_bar_h, yc_el_h + y_bar_h + 1] fl_1 = [fl_el , fl_el , fl_el ] xc_2 = [xc_el_v + x_bar_v - 1, xc_el_v + x_bar_v, xc_el_v + x_bar_v + 1] yc_2 = [yc_el_v + y_bar_v - 5, yc_el_v + y_bar_v, yc_el_v + y_bar_v + 5] fl_2 = [fl_el , fl_el , fl_el ] xc_barre = [xc_1, xc_2] yc_barre = [yc_1, yc_2] fl_barre = [fl_1, fl_2] * 0.15 ; * 0.15 car nbreuses sources. ; Sources larges (7 * 7 pixels) x_el = [-3, -2, -1, 0, 1, 2, 3, -3, -2, -1, 0, 1, 2, 3, $ -3, -2, -1, 0, 1, 2, 3, -3, -2, -1, 0, 1, 2, 3, $ -3, -2, -1, 0, 1, 2, 3, -3, -2, -1, 0, 1, 2, 3, $ -3, -2, -1, 0, 1, 2, 3] y_el = [-3, -3, -3, -3, -3, -3, -3, -2, -2, -2, -2, -2, -2, -2, $ -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, $ 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, $ 3, 3, 3, 3, 3, 3, 3] f_el = [ 0, 0, 3, 3, 3, 0, 0, 0, 3, 3, 3, 3, 3, 0, $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, $ 3, 3, 3, 3, 3, 3, 3, 0, 3, 3, 3, 3, 3, 0, $ 0, 0, 3, 3, 1, 0, 0] xc_lar = [xc_lar_1 + x_el, xc_lar_2 + x_el] yc_lar = [yc_lar_1 + y_el, yc_lar_2 + y_el] fl_lar = [f_el , f_el ] * 0.5 ; sources nbreuses. ; Sources eparses xc_ep = [40, 43, 46, 49, 52, 55, 60, 64, 67, 70, 73, 75, 78, 81, 84] yc_ep = [55, 65, 48, 54, 75, 64, 48, 75, 45, 67, 78, 50, 66, 73, 59] fl_ep = [ 4, 8, 0, 3, -2, 5, -2, 7, 5, -3, 4, 6, -3, 6, 4] * 1 if (i_simple eq 1) then begin a = 2 ; echelle caracteristique (pixels). xc_ep = [-a, -a, 0, a, a, 2*a] + 64 yc_ep = [-a, a, 0, -a, a, 2*a] + 64 fl_ep = [ 3, 3, 3, 3, 3, 0] endif ; Source centrale tres intense xc_int = [-1, 0, 1, -1, 0, 1, -1, 0, 1] + 64 ;48 yc_int = [-1, -1, -1, 0, 0, 0, 1, 1, 1] + 64 ;48 fl_int = [ 0, 0, 0, 0, 16, 0, 0, 0, 0] ; fl_int = [ 0, 0, 0, 0, 16, 0, 0, 0, 0] ; source ponctuelle ; fl_int = [ 2, 2, 2, 2, 4, 2, 2, 2, 2] ; coins non a zero. ; fl_int = [ 0, 2, 0, 2, 4, 2, 0, 2, 0] ; coins a zero. ; Pour une source tres compacte : ; 2 ; 242 ; 2 ; Usuel : xc_int=64, yc_int=60 ; fl_int= [ 2, 2, 2, 2, 2, 2, 2, 2, 2] ; (on peut mettre les coins a zero). ; FIN DES PARAMETRES FIXES (OU PLUS RAREMENT REGLABLES) DU SOLEIL ; Ensemble des sources plus ou moins compactes (flilament + barres + etc.) xc = [xc_fil , xc_barre , xc_ep, xc_lar, xc_int] yc = [yc_fil , yc_barre , yc_ep, yc_lar, yc_int] f_c = [amul_fil*fl_fil, amul_barre*fl_barre*2, amul_ep*fl_ep, $ amul_lar*fl_lar, amul_int*fl_int] ; CALCUL DES DISTRIBUTIONS DE BRILLANCE ; Calcul de la distribution de brillance de la source large. ; On utilise une repartition du type : ; amp = (1 - r^p)^p pour r<1 (1) ; amp = 0 pour r>1 ; ou r est un rayon vecteur normalise (max = 1) et p un exposant. ; Dans une distribution anisotrope elliptique les lignes de niveau ; sont des ellipses homothetiques d'equation ; x^2/a^2 + y^2/b^2 = 1 (2) ; On pose a = m * A et b = m * B ou A et B sont les 1/2 axes de ; l'ellipse la plus grande (de niveau nul). On a donc ; x^2/A^2 + y^2/B^2 = m^2 (3) ; Le 1er membre de (3) est donc constant sur une ellipse donnee et on ; peut le prendre comme "rayon vecteur" r dans l'eq (1) ; Definition du champ reduit carre dans lequel la souce large est centree ; La source etant large par hypothese, on la definit directement sur ; np2*np2 pts, sans interpolation. On la calcule sur un champ reduit, ; centre au centre su soleil calme, et dont la largeur est limitee ; par le bord du champ le plus proche. np2 = 2 * min([npi-xb, xb, npi-yb, yb]) ; Calcul des coordonnees centrees et reduites, x/A et y/B et du "rayon vec- ; teur r limit'e a [0, 1]. x_r = (findgen(np2) - np2/2) / larg_ew ; x/A notation precedente. y_r = (findgen(np2) - np2/2) / larg_ns ; y/B ------------------- mx = x_r # replicate(1, np2) my = replicate(1, np2) # y_r re = sqrt(mx^2 + my^2) ; Definitions domaines interieurs et exterieurs a la bordure du Soleil calme dom_ext = where(re gt 1, ic_ext) ; domaine exterieur a la limite du SC. dom_int = where(re le 1, ic_int) ; ------- interieur ----------------- ; Limitations a de x_r, t_r et re if (ic_ext gt 0) then re(dom_ext) = 1 amp = (1 - re^ae)^ae ; Immersion dans un tableau npi*npi . ima_b = fltarr(npi, npi) ; b comme base. ima_b(npi/2-np2/2, npi/2-np2/2) = flux_b / total(amp) * amp ; Definitions utiles pour la contribution genre CME max_SC = max(ima_b) dom_05 = where(ima_b gt 0.5*max_SC) ; domaine interieur "notable". ; Controle de source large simulee au centre du champ total (npi/2, npi/2). icon = 0 if (icon eq 1) then begin window, 0 & wset, 0 im_visu = congrid(ima_b, 512, 512, cubic=-0.5) shade_surf, im_visu xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'source large simulee au centre du champ total' window, 1 & wset, 1 tvscl, im_visu stop endif ; Images de la source large a sa position finale. im_b = shift(ima_b, xb - npi/2, yb - npi/2) ; b comme base. ; Controle de source large simulee a sa position finale. icon = 0 if (icon eq 1) then begin help, im_b ; Traces 2D window, 0 & wset, 0 shade_surf, im_b xyouts, 0.5, 0.07, /normal, alignment=0.5, $ 'Base large sur tout le champ' window, 1 & wset, 1 tvscl, congrid(im_b, 512, 512, cubic=-0.5) xyouts, 0.5, 0.07, /normal, alignment=0.5, $ 'Base large sur tout le champ' ; Traces 1D absc = - npi/2 + indgen(npi+1) window, 2 & wset, 2 plot, absc, im_b(* , npi/2), xstyle=1 oplot, absc, im_b(npi/2, * ), linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Source large position finale : coupe EW (___) et NS(...)' stop endif ; Calcul de la distribution de brillance des sources compactes ; Se fait en plusieurs etapes : ; - calcul d'une composante situee au centre d'un champ reduit. Toutes ; les sources elementaires etant centrees sur des pixels entiers, ; il n'est pas utile de faire une forte interpolation comme dans ; SIMUL_VISIB_1 ou on veut pouvir regler precisement les positions ; des 2 composantes du pseudo-Cygne. ; - calcul des 2 composantes de la source compacte avec la meme inter- ; polation ; - reduction son centre ; de gravite au milieu d'un champ carre comprenant np1 pixels de part ; et d'autre de son centre, c'est a dire s'etendant sur une grille de ; (2*np1 + 1) * (2*np1 + 1) pixels. On prends le centre de ce champ ; reduit en (npi/2, npi/2) et on choisit np1 assez petit pour eviter ; les underflows dans le calcul des exponentielles. ; - pair (precaution du 9 oct 00 pour eviter des decalages intempes- ; tifs lors de l'immersion) ; On decrit la source compacte par 128*128 pts sur ce champ reduit (in- ; terpolation par un facteur 128/np1). ; Calcul d'une source compacte dans un champ reduit interpol'e np1 = fix(3.0*larg_c + 0.5) ; 3 *larg_c pour ne pas couper ; les pieds des gaussiennes. if (np1 gt npi/4) then np1 = npi/4 f_i = 1 ; 16 ; facteur d'interpolation (vieillerie). np3 = np1 * f_i ; 2*np3 + 1 nbre de points sur le ; champ reduit interpole. d_c = shift(dist(2*np3 + 1), np3, np3) / f_i ; d_c distance au centre dans le champ reduit interpole, expri- ; mee en canaux initiaux. sig = larg_c / sqrt(2 * alog(2)) ; ecart type de la gaussienne ; exp(-x^2 / (2 * sig^2)) arg = - d_c^2 / (2 * sig^2) ; argument gaussienne. amp = exp(arg) ; composante elementaire, decrite par ; (2*np1 + 1) * (2*np1 + 1) pts sur ; le champ reduit. ; Controle du bon positionnement de amp icon = 0 if (icon eq 1) then begin absc = (- np3 + indgen(2 * np3 + 1)) / float(f_i) ; en canaux initiaux. window, 0 & wset, 0 plot, absc, amp(*, np3) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Coupe EW d une composante de source compacte' window, 1 & wset, 1 plot, absc, amp(np3, *) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Coupe NS d une composante de source compacte' stop ; Conclusion (6 oct 2000) : le max est bien en 0, cad au milieu. endif ; Immersion de amp dans un champ complet interpol'e npint = npi * f_i amp_i = fltarr(npint, npint) amp_i(npint/2 - np3, npint/2 - np3) = amp ; Controle de amp_i icon = 0 if (icon eq 1) then begin toto = amp_i(npint/2-np3:npint/2+np3, npint/2-np3:npint/2+np3) absc = (- np3 + indgen(2 * np3 + 1)) / float(f_i) ; en canaux initiaux. window, 0 & wset, 0 plot, absc, toto(*, np3), xstyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Coupe EW source elementaire immergee (partie centrale) window, 1 & wset, 1 plot, absc, toto(np3, *), xstyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Coupe NS source elementaire immergee (partie centrale) stop endif ; Reduction a npi*npi pts d'une source compacte unite en milieu de champ amp_npi = congrid(amp_i, npi, npi, cubic=-0.5) ; Le facteur de reduction etant entier, congrid prend 1 pt sur f_i pts. amp_i = 0 ; destruction de amp_i. ; Calcul de la brillance de l'ensemble des sources compactes im_c = fltarr(npi, npi) imax = n_elements(xc) for i=0, imax-1 do begin im_el = shift(amp_npi, xc(i)-npi/2, yc(i)-npi/2) im_c = im_c + im_el * f_c(i) / total(im_el) endfor icon = 0 if (icon eq 1) then begin window, 0 im_c_512 = congrid(im_c, 512, 512, cubic=-0.5) tvscl, im_c_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Toutes les sources compactes' window, 1 shade_surf, im_c_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Toutes les sources compactes' stop endif ; Calcul distribution de brillance arche genre CME. x_CME = x_CME * larg_ew y_CME = y_CME * larg_ns ar_0 = shift(dist(npi), npi/2, npi/2) ; distance au centre du champ. ar_0(where(ar_0 eq 0)) = 0.5 ; eviter division par zero plus bas. ; Calcul du domaine annulaire ar = shift(ar_0, x_CME, y_CME) ; dist radiale a partir centre arc CME ; ar = ar_0 ; ---------------------------- disque dom_e = where(ar le r_CME) dom_i = where(ar le r_CME - dr_CME) tab_e = intarr(npi, npi) tab_i = intarr(npi, npi) tab_e(dom_e) = 1 tab_i(dom_i) = 1 dom_a = where( (tab_e eq 1) and (tab_i eq 0) ) ; Anneau CME uniforme im_CME = fltarr(npi, npi) im_CME(dom_a) = intens_CME * max_SC ; Determination du domaine ou le CME s'enracinera a la base du soleil calme val_rel = 0.5 ; valeur relative souhaitee de la base ; a la racine du CME. val = val_rel * max_SC ; ------ absolue --------- ... f_e = 1 ; facteur d'elargissement du domaine ; de valeurs sur le flanc du SC. n_e = 1 ; pour controle a1: dom_rac = where( (im_b le 1.01*f_e*val) and (im_b ge 0.99/f_e*val) and $ (im_CME ne 0), ic) ; domaine souhaite pour "l'enracine- ; ment" du CME sur le flanc du SC. if (ic lt 2) then begin ; il faut au moins 2 elements pour la ; fonction "moment" f_e = 1.01 * f_e ; elargissement progressif si le domai- ; ne est vide. n_e = n_e + 1 goto, a1 endif mom_ar_rac = moment(ar_0(dom_rac)) ; ensemble des distances au centre de ; la region d'enracinement. ar_rac = mom_ar_rac(0) ; moyenne du precedent. ; Ponderation du CME par une fonction radiale decroissante a partir du centre ; du champ, nulle sur le disque et valant 1 a la distance de l'enraci- ; nement pond = (ar_rac / ar_0) ^ puiss ; Mise a zero de la ponderation du CME aux distances de centre du champ infe- ; rieures a celle de l'enracinement dom_0 = where(pond gt 1) pond(dom_0) = 0 im_CME = im_CME * pond ; Raccord continu du CME au Soleil calme (on manipule des brillance en unites ; du maximum de la base max_SC) dom = where(im_CME ne 0) ; pour eviter division par zero ci-apres. im_CME_npi = im_CME ; sauvegarde non ponderee. fac_rac = (im_CME(dom)/intens_CME - im_b(dom)/val) / $ (im_CME(dom)/intens_CME) ; fac_rac varie de 0 a 1 sur le flanc de la base. dom_neg = where((fac_rac lt 0), ic) if(ic gt 0) then fac_rac(dom_neg) = 0 ; pour eviter que le CME soit < 0 pres de son enracinement. im_CME(dom) = amul_CME * im_CME(dom) * fac_rac icon = 0 if (icon eq 1) then begin ; toto = fltarr(npi, npi) ; toto(dom_0) = 1 ; window, 0 ; shade_surf, toto ; xyouts, 0.5, 0.97, /normal, alignment=-0.5, $ ; 'Domaine de mise a zero de pond' ; window, 1 ; shade_surf, pond ; xyouts, 0.5, 0.97, /normal, alignment=-0.5, $ ; 'Fonction de ponderation de CME' ; stop im_CME_npi_512 = congrid(im_CME_npi, 512, 512, cubic=-0.5) window, 0 tvscl, im_CME_npi_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'CME non pondere' window, 1 shade_surf, im_CME_npi_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'CME non pondere' im_CME_512 = congrid(im_CME, 512, 512, cubic=-0.5) window, 2 tvscl, im_CME_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'Source genre arche CME' window, 3 shade_surf, im_CME_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'Source genre arche CME' stop endif ; Source finale (base large + filament + composantes compactes + CME) : ; On interdit a l'objet final d'etre negatif (ce qui peut arriver du ; fait du couloir de filament si amul_comp est grand. bmul = 0.01 ; 0.01 pour reduire le flux a une valeur ; realiste de SFU. soleil = bmul * (im_b + im_c + im_CME) dom_neg = where(soleil lt 0, ic) if (ic gt 0) then begin print, 'ATTENTION objet <0 sur filament, remis a zero' soleil (dom_neg) = 0 endif icon = 0 if (icon eq 1) then begin soleil_512 = congrid(soleil, 512, 512, cubic=-0.5) window, 0 & wset, 0 tvscl, soleil_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Soleil calme simule' window, 1 & wset, 1 shade_surf, soleil_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Soleil calme simule' stop endif flux_final = total(soleil) ; Impression des flux des composantes et du flux total de l'objet ; le flux de la base flux_b est deja calcul'e. flux_b = bmul * flux_b flux_compact = bmul * total(im_c) flux_CME = bmul * total(im_CME) flux_total = flux_b + flux_compact + flux_CME print, ' ' print, 'SIMUL_VISIB_2 (soleil complexe) : ' + $ 'flux des objets composants :' print, ' flux de la base =', flux_b print, ' com compactes =', flux_compact print, ' CME =', flux_CME print, ' total objet =', flux_total print, ' ' ; CALCUL VISIBILITE "QUASI_REELLE" CENTREE tf_obj = fft(shift(soleil, -npi/2, -npi/2), 1) ; "quasi-reelle. visib = shift(tf_obj, npi/2, npi/2) ; centree dans espace TF. visib(npi/2, npi/2) = 0 ; simuler comp continue obser- ; vee n'est pas acquise. ; Controle icon = 0 if (icon eq 1) then begin print, 'SIMUL_VISIB_2 : visualisation soleil simule' soleil_512 = congrid(soleil, 512, 512, cubic=-0.5) window, 4 shade_surf, soleil_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'SIMUL_VISIB_2 : Soleil reel simule' write_gif, 'Sol_simule.gif', bytscl(soleil_512) window, 4 tvscl, soleil_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'SIMUL_VISIB_2 : Soleil reel simule' stop endif end ; fin de SIMUL_VISIB_2 ;______________________________________________________________________________ ; SIMUL_HARM_N est appele par RH_DPATCHFITS_NRH ou par RH_SYNTHESE, toujours ; pour une simulation, soit pour : ; . Soleil complexe (non polar ou polar), si i_choix_sim = 2 . ; Dans ce cas i_calib (0 ou 3) est indifferent et ne figure pas dans ; l'appel a SIMUL_VISIB_2 . ; . Cygne a deux composantes pour la calibration si i_choix_sim = 1 . ; Dans ce cas i_calib = 1 ou 2 et devient direct_auto pour l'appel a ; SIMUL_VISIB_1 a travers VISIB_CYG (calibration directe) ou directement ; par SIMUL_HARM_N (calibration auto). PRO SIMUL_HARM_N, $ ; entree i_choix_sim, i_calib, freq, t_moy, correl, $ i_pave_plein, i_sys_lin, n_bord, $ ; controles seuls. ; sortie gain, phase, $ ; gain et phases des antennes objet, pix_obj, $ ; objet et pixel sur l'objet visib_rh, mat_u, mat_v, $ ; visib sur la grille RH, coord des ; pts de la grille. harm_N ; vecteur simule des donnees de Nancay. ; visib : TF (centree dans espace TF) de l'objet simul'e en sa position reelle. ; But: produire un vecteur de 576 correlations harm_N correspondant a une sour- ; ce de structure connue observee avec un reseau ayant des gains et phases ; connus. On simule harm_N pour des dates posterieures au 8 oct 98. ; Appel : DPATCHFITS, RH_SYNTHESE ; Rem: "attempt to call an undefined procedure" quand le nom du ficher etait ; simul_harm_N.pro : lors de l'appel, unix convertissait le nom de fichier ; en minuscules et ne trouvait pas simul_harm_N. ; => ne pas melanger minuscules et majuscules dans les noms de fichiers. ; Modifications: ; 99 dec 01 debut ecriture ; 00 avr 06 ajout d'une base large a l'image simulee ; 00 avr 17 la position de l'ext NS est changee en 44. NS46 devient NS44. ; -- --- 25 -------------------------------------- 45 NS44 ------- NS45. ; mai 2 verification systematique des modifications subsequentes. ; 01 sep 13 . la simulation de l'objet soleil est reportee dans les sous- ; procedures SIMUL_VISIB_1 et SIMUL_VISIB_2, en tete de ; fichier. ; . introduction de i_choix_sim dans la liste d'entree pour choi- ; sir le type de source simulee : ; . Cygne pour une calibration (directe ou auto) ; . Soleil complexe pour des essais de CALC_IMAGE_HELIO (test ; de clean etc.) et pour CALIBRATION. ; dans les 2 cas on peut simuler un reseau dephase. ; . n_bord (hauteur des bandes N et S pour ajuster base dans ; RH_MALC_IM_2D) est defini dans DPATCHFITS et passe a toutes ; les procedures appelantes de RH_MALC_IM_2D. ; 21 ajout de l'objet simule dans la liste en sortie, a des fins de ; controle dans WCLEAN. ; 25 ajout de i_pave_plein a la liste en entree, pour passage a ; RH_MALC_IM_2D pour controles seulement. ; 02 fev 16 "gain" et "phase" en sortie pour comparaison avec amul et bmul ; dans CALIBRATION. ; 22 visib (TF de l'objet simul'e) en sortie pour controles dans ; CALIBRATION. ; mar 28 ajout de i_calib freq et t_moy dans la liste d'appel pour ; pouvoir simuler des donnees pour une calibration directe : ; les donnees etant comparees a la visibilte du modele ana- ; lytique simple du Cygne a une heure moyenne entre le debut ; et la fin choisis dans RH_DPNEW, il faut qu'elles soient ; produites a partir du meme modele. ; 03 sep 3 usage de i_pave_plein=2 pour signifier, a partir de VIS_MOD_ ; CYGNE, qu'il faut : ; . remettre i_pave_plein a zero, ; . transmettre i_champ_uv=2 a VISIB_CYG, qui le transmet ; a SIMUL_VISIB_1. ; Ceci evite de retoucher la liste d'entree de SIMUL_HARM_N, ; qui est appel'e de plusieurs programmes. ; 03 oct 10 extension au cas des antennes anti-aliasing AA. On reconnait ; ce cas grace a la dimension de correl. ; 04 nov 9 mise au clair des comMentaires sur SIMUL_VISIB_1, VISIB_CYG, ; ajout de pix_obj, mat_u, mat_v a la liste de sortie. ; Notations ; correl tableau (2, 576) indiquant les numeros des 2 antennes ; correspondant a chacun des 576 harmoniques. Permet de faire ; le passage de harm_N a harm. ; harm_N complexarr (576) harmoniques, rang'es au format Nancay. Rempli ; en non polar ou polar juste avant l'appel a calc_image_helio. ; harm complexarr (128, 128) tableau des harm dans le plan (u, v). ; i_choix_sim 1 Cygne, 2 Soleil ; i_calib type de calibration : directe avec modele donn'e (1) ou auto ; avec modele resultant du lissage des donnees (2). ; objet objet simul'e en sa position reelle (taille 512*512 pour le ; Cygne avec i_calib=1). ; visib_rh TF (centree en milieu de tableau dans l'espace TF) de l'objet ; simul'e (pour lequel l'origine des coordonnees est en milieu ; de tableau) en sa position reelle. ; . i_calib = 1 (calibration directe: ; - objet defini dans l'espace reeel, et pix_obj en rad^-1 ; taille pour le Cygne 512*512 ; - visib_rh est la visibilite sur la grille (completee) ; des harmoniques du RH (taille 128*128). ; - mat_u et mat_v sont les tableaux des coord uv de ces ; points, en rad^-1. ; . i_calib = 2 : ; - "objet" est defini sur une grille a maille carree dans ; l'espace interferometrique, ; - visib_rh est definie sur la grille carre conjuguee, ; - mat_u et mat_v produiotes simplement avec "indgen" ; et "replicate". ; Rangement des harmoniques dans harm_N ; NS0 avec E2, E1, H1 a H16 0 a 17 ; NS1 ---- ------------------------------ 18 a 35 ; NS8 ---- ------------------------------ 144 a 161 ; NS23 ---- ------------------------------ 414 a 431 ; NS45 ---- ------------------------------ 432 a 449 ; jusqu'au 7 oct 98 inclus: ; E0 ---- NS0 a NS17 450 a 467 ; -- ---- NS10 a NS23 468 a 481 ; -- ---- NS45 482 ; -- ---- NS7 a NS9 483 a 485 ; a partir du 8 oct 98 inclus: ; NS45 ---- NS0 a NS17 450 a 467 ; -- ---- NS10 a NS23 468 a 481 ; -- ---- NS45 482 ; -- ---- NS7 a NS9 483 a 485 ; E2 ---- E2, E1, H1 a H16 486 a 503 ; E1 ---- E2, E1, H1 a H16 504 a 521 ; E0 ---- E2, E1, H1 a H16 522 a 539 ; NS8 ---- NS0 a NS17 540 a 557 ; NS0 ---- NS10 a NS23 et NS45 558 a 572 ; --- ---- NS7 a NS9 573 a 575 ; AA1 ---- NS0 a NS17 576 a 593 ; AA2 ---- ---------- 594 a 611 ; AA3 ---- ---------- 612 a 629 ; AA4 ---- ---------- 630 a 647 ; XX1 ---- YY1 648 a 665 ; XX2 ---- YY2 666 a 683 ; XX3 ---- YY3 684 a 701 ; XX4 ---- YY4 702 a 719 ; Rangement des antennes (comme dans les tableaux des coord des antennes) ; num ant bew num ant bew num ant num ant ; 0 E2 -48 3 H1 -14 19 NS 0 35 NS 16 ; 1 E1 -16 4 H2 -12 20 NS 1 36 NS 17 ; 2 E0 -15 5 H3 -10 21 NS 2 37 NS 18 ; 6 H4 -8 22 NS 3 38 NS 19 ; 7 H5 -6 23 NS 4 39 NS 20 ; 8 H6 -4 24 NS 5 40 NS 21 ; 9 H7 -2 25 NS 6 41 NS 22 ; 10 NS0 0 26 NS 7 42 NS 23 ; 11 H9 2 27 NS 8 43 NS 45 ; 12 H10 4 28 NS 9 ; 13 H11 6 29 NS 10 44 AA 1 ; 14 H12 8 30 NS 11 45 AA 2 ; 15 H13 10 31 NS 12 46 AA 3 ; 16 H14 12 32 NS 13 47 AA 4 ; 17 H15 14 33 NS 14 ; 18 H16 16 34 NS 15 ; Discernement des cas sans ou avec les antennes AA taille = size(correl) nbre_harm = taille(2) ; taille = [nbre de dim, dim1, dim2, ...] if (nbre_harm eq 576) then nbre_ant = 44 if (nbre_harm gt 576) then nbre_ant = 48 ; Gestion de la rustine sur i_pave_plein pour l'usage avec VIS_MOD_CYGNE (dans ; fichier du meme nom). i_champ_uv = 1 if (i_pave_plein eq 2) then begin i_pave_plein = 0 ; remise en ordre pour RH_MALC_IN_2D. i_champ_uv = 2 ; pour usage special a partir de VIS_MOD_CYGNE ; (dans fichier du meme nom) pour aller assez ; loin pour ajuster la largeur larg_c des ; sources du modele. endif if (i_pave_plein eq 3) then begin i_pave_plein = 0 ; remise en ordre pour RH_MALC_IN_2D. i_champ_uv = 1 ; quand on n'a pas besoin de decrire le 2eme ; max. endif np = 128 ; dimension de l'objet simule (identique au npi de dpnew). ; Definition des facteurs de gain compelexe des antennes. Utilises seulement en ; fin de procedure, ils figurent au debut par commodite. gain = fltarr(nbre_ant) ; tableau des gains reels. phase = fltarr(nbre_ant) ; tableau des phases. gain_c = complexarr(nbre_ant) ; gains complexes produits avec gains ; et phases. gain_unif = [ $ 1.00, 1.00, 1.00, $ ; Ext 2 1 0 1.01, 1.01, 0.99, 0.99, 1.01, 1.01, 0.99, $ ; H 1 a H7 1.01, $ ; NS0_ew H8 0.99, 1.01, 0.99, 1.01, 0.99, 1.01, 0.99, 1.01, $ ; H 9 a 16 0.99, 0.99, 1.01, 1.01, 0.99, 0.99, 1.01, 1.01, 0.99, 0.99, $ ; NS 0 a 9. 1.00, 1.01, 0.99, 0.99, 1.01, 1.01, 0.99, 0.99, 1.01, 1.01, $ ; NS 10 a19 1.00, 1.00, 1.00, 1.00, 1.20, $ ; NS 20 23, ; NS 45. 1.00, 1.00, 1.00, 1.00] ; AA1 a 4. phase_unif = [ $ ; en degres. 00, -00, -0, $ ; Ext 2 1 0 01, 01, -01, -01, 01, 01, -01, $ ; H 1 a 7 -0, $ ; NS0_ew H8 -01, 01, 01, -01, -01, 01, 01, -01, $ ; H 9 a 16 -00, 01, -01, -01, 01, 01, -01, -01, 01, 01, $ ; NS 0 a 9. -01, -01, 01, 01, -01, -01, 01, 01, -01, -01, $ ; NS 10 a 19 -01, -01, 01, 01, -20, $ ; NS 20 a 23 et NS 45 -01, -01, 01, 01 ] ; AA1 a AA4 gain_regulier = [ $ 1.10, 0.90, 1.10, $ ; Ext 2 1 0 0.90, 1.10, 0.90, 1.10, 0.90, 1.10, 0.90, $ ; H 1 a H7 1.10, $ ; NS0_ew H8 0.90, 1.10, 0.90, 1.10, 0.90, 1.10, 0.90, 1.10, $ ; H 9 a 16 1.20, 0.90, 1.10, 0.90, 1.10, 0.90, 1.10, 0.90, 1.10, 0.90, $ ; NS 0 a 9. 1.10, 0.90, 1.10, 0.90, 1.10, 0.90, 1.10, 0.90, 1.10, 0.90, $ ; NS 10 a19 1.10, 0.90, 1.10, 0.90, 1.10, $ ; NS 20 23, ; NS 45. 1.10, 0.90, 1.10, 0.90 ] ; AA1 a 4. phase_regulier = [ $ ; en degres. 10, -10, 10, $ ; Ext 2 1 0 10, -10, 10, -10, 10, -10, 10, $ ; H 1 a 7 -10, $ ; NS0_ew (H8) -10, 10, -10, 10, -10, 10, -10, 10, $ ; H 9 a 16 ; 10, -10, 10, -10, 10, -10, 10, -10, $ ; H 9 a 16 -20, 10, -10, 10, -10, 10, -10, 10, -10, 10, $ ; NS0_ns a NS9. -10, 10, -10, 10, -10, 10, -10, 10, -10, 10, $ ; NS 10 a 19 -10, 10, -10, 10, 20, $ ; NS 20 a 23 et NS 45 -10, 10, -10, 10 ] ; AA1 a AA4 gain_aleatoire = [ $ 1.10, 1.10, 0.80, $ ; Ext 2 1 0 1.07, 0.95, 0.92, 1.12, 0.93, 1.06, 1.03, $ ; H 1 a H7 1.13, $ ; NS0_ew H8 1.12, 0.93, 0.92, 1.07, 1.09, 0.84, 1.12, 0.95, $ ; H 9 a 16 0.85, 1.07, 0.90, 1.05, 0.94, 0.85, 1.13, 1.13, 1.17, 0.80, $ ; NS 0 a 9. 1.05, 0.95, 1.07, 1.03, 0.93, 0.86, 1.05, 1.04, 0.95, 1.04, $ ; NS 10 a19 0.98, 1.12, 0.90, 0.91, 1.05, $ ; NS 20 23. ; NS 45. 0.98, 1.12, 0.90, 0.91 ] ; AA1 a 4. phase_aleatoire = [ $ ; en degres. 10, -8, 7, $ ; Ext 2 1 0 -10, 7, -5, 12, 17, -10, 8, $ ; H 1 a 7 7, $ ; NS0_ew H8 -12, -13, -10, 11, -4, 3, -14, -2, $ ; H 9 a 16 00, 13, -12, 15, 10, 6, 12, -15, -12, 8, $ ; NS 0 a 9. -9, 12, 11, -15, -14, 20, 13, -16, -10, 13, $ ; NS 10 a 19 2, -5, 14, -10, 15, $ ; NS 20 a 23 et NS 45 2, -5, 14, -10 ] ; AA1 a AA4 ; phase_aleatoire = [ $ ; en degres. ; -10, -10, 07, $ ; Ext 2 1 0 ; -15, 10, -5, -0, 12, -8, 17, $ ; H 1 a 7 ; 14, $ ; NS0_ew H8 ; -12, -10, -15, 12, -17, 10, -18, -15, $ ; H 9 a 16 ; 0, 18, -11, 25, 8, 14, 45, -26, -27, 16, $ ; NS 0 a 9. ; -19, 21, 28, -31, -34, 18, 23, -27, -66, 19, $ ; NS 10 a 19 ; 28, -15, 34, -13, 25, $ ; NS 20 23 et NS 45 ; 28, -15, 34, -13 ] ; AA1 a AA4. ; Choix du type de gain (aleatoire ou regulier) et reduction des inegalites de ; gain et phase pour mise au point. choix_gain = 1 ; 0 uniforme, 1 regulier, 2 aleatoire. amul_gain = 0.001 ; 4.9 limite comportement aberrant. amul_phi = 0.001 ; Calcul de la visibilite simulee if (i_choix_sim eq 1) then begin ; simulation Cygne pour une calibra- ; tion, appel de SIMUL_VISIB_1, soit ; . directement (calibration auto), ; . a travers VISIB_CYG (calibration ; directe) direct_auto = i_calib ; pour concordance de notation if (i_calib eq 1) then begin ; calcul de visib sur grille reelle ; (non carree) et complete du RH pour ; une calibration directe. VISIB_CYG, $ ; procedure dans CALIB_SUB. ; VISIB_CYG fixe lui meme direct_auto = 1 et ; appelle SIMUL_VISIB_1 ; entree i_champ_uv, $ freq, t_moy, $ ; freq en MHz correl, $ ; ne sert qu'a un calcul supprimable. ; sortie objet, $ ; objet simule (transmis de SIMUL_VISIB_1), ; dans l'espace reel (cosddh, ddelta) pour ; une calibration directe. Il est defini ; par 512 points sur un champ de 22' d'arc ; dans SIMUL_VISIB_1 . Ne sert qu'a des ; controles dans SIMUL_HARM_N. pix_obj, $ ; pixel sur l'objet, (rad) visib_rh, $ ; visibilite (centree en milieu de tableau) ; sur la grille completee des harm du RH, ; c'est a dire une grille 128*128, avec une ; maille plus ou moins en parallellogramme. mat_u, mat_v ; matrices des coordonnees des points ; (rad^-1) ou est donnee visib_RH. endif if (i_calib eq 2) then begin ; calibration de type auto. direct_auto = i_calib ; rappel mnemotechnique. i_champ_uv = 1 ; bidon mais laisser a 1 SIMUL_VISIB_1, $ ; entree direct_auto, $ ; sortie ; np, $ ; dimension image. objet, $ ; objet simul'e en position relle dans ; l'espace interferometrique pour ; une calibration auto pix_obj, $ ; utilise seulement quand SIMUL_VISIB_1 ; est appele par VISIB_CYG. Ici inu- ; tile (=1) visib_rh, $ ; TF (centree en milieu de tableau) de ; l'objet simul'e en position reelle. pix_uv ; pixel sur visib_rh. ici inutile (=1). taille = size(visib_rh) np = taille(1) visib_rh(np/2, np/2) = 0 ; pour simuler le fait que la composan- ; te continue n'est pas acquise par ; le RH. vec_u = indgen(128) - 64 v_rep = replicate(1, 128) mat_u = vec_u # v_rep ; mat_u et mat_v pour compatibilite mat_v = v_rep # vec_u ; avec le cax i_calib=1 . endif endif if (i_choix_sim eq 2) then begin ; simulation pour CALC_IMAGE_HELIO. ; print, 'SIMUL_HARM_N : i_choix_sim =', i_choix_sim ; print, '------------ : simulation pour CALC_IMAGE_HELIO' SIMUL_VISIB_2, $ ; definitions source dans SIMUL_VISIB_2 ; sortie np, $ ; dimension image. objet, $ ; objet simul'e en position relle. visib_rh ; TF (centree dans espace TF) d'objet simul'e ; en position reelle. endif ; Remplissage du vecteur harm_N npuv = i_champ_uv * np ; rappel : i_champ_uv=1 sauf dans l'utilisation avec VIS_MOD_CYGNE ; (dans fichier du meme nom), qui visualise la TF du modele de Cygne ; ou il vaut 2. harm_N = complexarr(nbre_harm) ; Rangement des harmoniques dans harm_N ; NS0 avec E2, E1, H1 a H16 0 a 17 ; NS1 ---- ---------------- 18 a 35 ; NS8 ---- ---------------- 144 a 161 ; NS23 ---- ---------------- 414 a 431 ; NS45 ---- ---------------- 432 a 449 ; jusqu'au 7 oct 98 inclus: ; E0 ---- NS0 a NS17 450 a 467 ; -- ---- NS10 a NS23 468 a 481 ; -- ---- NS45 482 ; -- ---- NS7 a NS9 483 a 485 ; a partir du 8 oct 98 inclus: ; NS45 ---- NS0 a NS17 450 a 467 ; -- ---- NS10 a NS23 468 a 481 ; -- ---- NS45 482 ; -- ---- NS7 a NS9 483 a 485 ; E2 ---- E2, E1, H1 a H16 486 a 503 ; E1 ---- E2, E1, H1 a H16 504 a 521 ; E0 ---- E2, E1, H1 a H16 522 a 539 ; NS8 ---- NS0 a NS17 540 a 557 ; NS0 ---- NS10 a NS23 et NS45 558 a 572 ; --- ---- NS7 a NS9 573 a 575 ; AA1 ---- NS0 a NS17 576 a 593 ; AA2 ---- ---------- 594 a 611 ; AA3 ---- ---------- 612 a 629 ; AA4 ---- ---------- 630 a 647 ; XX1 ---- YY1 648 a 665 ; XX2 ---- YY2 666 a 683 ; XX3 ---- YY3 684 a 701 ; XX4 ---- YY4 702 a 719 ; Harmoniques de E2 avec les NS for i = 0, 23 do begin ; harm_N(10*i) = visib_rh(npuv/2-48, i) ; NS0 a NS23 -> 15 fev 00. ; harm_N(18*i) = visib_rh(npuv/2-48, i) ; NS0 a NS23 apres 15 fev 00. harm_N(18*i) = visib_rh(npuv/2-48, npuv/2 + i) ; NS0 a NS23 apres ; 18 avr 02. endfor harm_N(432) = visib_rh(npuv/2 - 48, 45) ; NS45 ; Harmoniques du "pave": les NS avec E1 a H18 ; NS0 avec E2, E1, H1 a H16 0 a 17 ; NS1 ---- ------------------------------ 18 a 35 ; NS8 ---- ------------------------------ 144 a 161 ; NS23 ---- ------------------------------ a 431 ; NS45 ---- ------------------------------ 432 a 449 for ins = 0, 23 do begin ; boucle NS0 a NS23. for iew = -16, 16, 2 do begin ; boucle EW, harmoniques pairs. harm_N(1 + 18*ins + iew/2 + 8) = visib_rh(npuv/2 + iew, npuv/2+ins) endfor endfor for iew = -16, 16, 2 do begin ; boucle EW sur NS45 harm_N(433 + iew/2 + 8) = visib_rh(npuv/2 + iew, npuv/2 + 45) endfor ; Harmoniques de E2 avec E2, E1, H1 a H19 A VOIR A VOIR A VOIR A VOIR ; E2 avec E2, E1, H1 a H16 486 a 503 harm_N(486) = visib_rh(npuv/2 , npuv/2) harm_N(487) = visib_rh(npuv/2 + 32, npuv/2) for iew = 34, 62, 2 do begin ; boucle de H1 a H15 (harm 64 interdit par Shannon). harm_N(488 + (iew-34)/2) = visib_rh(npuv/2 + iew, npuv/2) endfor ; Erreur avant le 4 mai 2000. Cette erreur produisait une mauvaise calibration ; de E2 avec H1 (les harmoniques dans harm_N etaient decales d'une ; unite, et comme la phase de la visibilite tourne vite cela produisait ; un decalage de phase d'environ 15 deg. ; for iew = 34, 63 do begin ; avant le 4 mai 2000 ; ; boucle de H1 a H15 (harm 64 interdit par Shannon). ; harm_N(488 + (iew-34)/2) = visib_rh(npuv/2 + iew, npuv/2) ; endfor ; Harmoniques de E1 avec E2, E1, H1 a H16. ; E1 avec E2, E1, H1 a H16 504 a 521 harm_N(504) = visib_rh(npuv/2 - 32, npuv/2) ; -16 jusqu'au 15 fev 00. harm_N(505) = visib_rh(npuv/2 , npuv/2) for iew = 2, 32, 2 do begin ; boucle de H1 a H16. harm_N(506 + (iew-2)/2) = visib_rh(npuv/2 + iew, npuv/2) endfor harm_1 = harm_N ; sauvegarde icon = 0 ; controle des harm de E1 if (icon eq 1) then begin window, 0 & wset, 0 hmax = max(abs(harm_N(504:521)), min=hmin) plot, yrange =[0, 1.2*hmax], ystyle=1, xstyle=1, $ title='abs(harm de E1 avec E2, E1, H1 a H16)', $ abs(harm_N(504:521)) stop endif ; Harmoniques de E0 avec E2, E1, H1 a H16 ; E0 avec E2, E1, H1 a H16 522 a 539 harm_N(522) = visib_rh(npuv/2 - 33, npuv/2) ; E2, ligne de base vers Est. harm_N(523) = visib_rh(npuv/2 - 1, npuv/2) ; E1, ligne de base vers Est. ; faute de signe (npuv/2 + 1) trouvee le 10 fev 2000. for iew = 1, 31, 2 do begin harm_N(524 + (iew-1)/2) = visib_rh(npuv/2 + iew, npuv/2) ; H1 a H16. endfor icon = 0 ; controle des harm de E0 if (icon eq 1) then begin window, 0 & wset, 0 absc = indgen(npuv) - npuv/2 ; tvscl, congrid(abs(visib_rh), 512, 512, cubic=-0.5) ; plot, absc, abs(visib_rh(npuv/2 - 34 : npuv/2 + 34, npuv/2)), $ title='harm EW purs de visib_rh' plot, absc, abs(visib_rh(*, npuv/2)), xrange=[-20, 20], xstyle=1, $ title='harm EW purs de visib_rh' window, 1 & wset, 1 plot, title='(harm_N(E0) sauf E2, __abs, - -R', $ abs (harm_N(523:539)) oplot, float (harm_N(523:539)), psym=2 ; 2 *, 3 ., 4 losange oplot, imaginary(harm_N(523:539)), psym=3 print, 'SIMUL_HARM_N : harm_N(523:527) (E0 avec E1 H1 H2 H3 H4' print, harm_N(523:527) stop endif icon = 0 ; juste pour visualiser la visibilite. if (icon eq 1) then begin window, 0 & wset, 0 tvscl, congrid(abs(visib_rh), 512, 512, cubic=-0.5) xyouts, /normal, 0.5, 0.97, 'abs(visib_rh)', alignment=0.5 window, 1 & wset, 1 diff = harm_N - harm_1 plot, diff(520:540), psym=2 ; 2 *, 3 ., 4 losange stop endif ; Harmoniques de NS8 avec NS0 a NS17 ; NS8 ---- NS0 a NS17 540 a 557 for nns = 0, 7 do begin ; nns num d'antenne NS. harm_N(540+nns) = visib_rh (npuv/2, npuv/2 + 8 - nns) ; NS0 a NS7. endfor harm_N(548) = visib_rh(npuv/2, npuv/2) ; NS8 avec NS8. for nns = 9, 17 do begin harm_N(549 + nns-9) = visib_rh (npuv/2, npuv/2-(nns-8)) ; NS9 a NS17. endfor ; Harmoniques de NS0 avec ; NS0 ---- NS10 a NS23 et NS45 558 a 572 ; --- ---- NS7 a NS9 573 a 575 for nns = 10, 23 do begin harm_N(558 + nns-10) = visib_rh(npuv/2, npuv/2 - nns) endfor harm_N(573) = visib_rh (npuv/2, npuv/2 - 7) harm_N(574) = visib_rh (npuv/2, npuv/2 - 8) harm_N(575) = visib_rh (npuv/2, npuv/2 - 9) ; Harmoniques avec NS45 ; NS45 ---- NS0 a NS17 450 a 467 ; -- ---- NS10 a NS23 468 a 481 ; -- ---- NS45 482 ; -- ---- NS7 a NS9 483 a 485 for nns = 0, 17 do begin harm_N(450 + nns) = visib_rh (npuv/2, npuv/2 + 45 - nns) endfor for nns = 10, 23 do begin harm_N(468 + nns - 10) = visib_rh (npuv/2, npuv/2 + 45 - nns) endfor harm_N(482) = visib_rh (npuv/2, npuv/2) for nns = 7, 9 do begin harm_N(483 + nns - 7) = visib_rh (npuv/2, npuv/2 + 45 - nns) endfor ; Harmoniques avec lea antennes AA if (nbre_harm gt 576) then begin for iaa = 1, 4 do begin visib_ns = visib_rh(npuv/2 - 1 - 2*(iaa -1 ), $ npuv/2 - 9 : npuv/2 + 9) ; coupe // axe v limitee a + ou - 9 visib_ns = reform (visib_ns) ; supprimer une dim. degeneree. visib_ns_rev = reverse(visib_ns) ; de v=+9 a v=-9. 19 elements harm_N(576 + 18 * (iaa - 1)) = visib_ns_rev(1:18) endfor endif ; fin du calcul des harmoniques des antennes AA. harm_N_theo = harm_N ; sauvegarde ; Fin de calcul de harm_N non perturbe ; Rangement des antennes (comme dans les tableaux des coord des antennes) ; num ant bew num ant bew num ant num ant ; 0 E2 -48 3 H1 -14 19 NS 0 35 NS 16 ; 1 E1 -16 4 H2 -12 20 NS 1 36 NS 17 ; 2 E0 -15 5 H3 -10 21 NS 2 37 NS 18 ; 6 H4 -8 22 NS 3 38 NS 19 ; 7 H5 -6 23 NS 4 39 NS 20 ; 8 H6 -4 24 NS 5 40 NS 21 ; 9 H7 -2 25 NS 6 41 NS 22 ; 10 NS0 0 26 NS 7 42 NS 23 ; 11 H9 2 27 NS 8 43 NS 45 ; 12 H10 4 28 NS 9 ; 13 H11 6 29 NS 10 44 AA 1 ; 14 H12 8 30 NS 11 45 AA 2 ; 15 H13 10 31 NS 12 46 AA 3 ; 16 H14 12 32 NS 13 47 AA 4 ; 17 H15 14 33 NS 14 ; 18 H16 16 34 NS 15 ; Choix des jeux de gains et phases complexes if (choix_gain eq 0) then begin gain = gain_unif & phase = phase_unif endif if (choix_gain eq 1) then begin gain = gain_regulier & phase = phase_regulier endif if (choix_gain eq 2) then begin gain = gain_aleatoire & phase = phase_aleatoire endif phase = amul_phi * phase gain = 1 + amul_gain * (gain - 1) ; Calcul du gain instrumental simul'e complexe phase = float(!pi/180 * phase) ; conversion radians. i = complex(0, 1) gain_c = gain * exp(i * phase) ; Controle du gain complexe icon = 0 if (icon eq 1) then begin n1 = 0 & n2 = 43 absc = n1 + indgen(n2 - n1) window, 0 & wset, 0 ymin = min(abs(gain_c(n1:n2)), max=ymax) dy = 0.05 * (ymax - ymin) g_n1n2 = gain_c(n1:n2) plot, absc, abs(g_n1n2), $ xstyle=1, yrange=[ymin-dy, ymax+dy], ystyle=1, $ title='module du gain' dom = where(float(g_n1n2) ne 0) arg = 180/!pi * atan(imaginary(g_n1n2), float(g_n1n2)) window, 1 & wset, 1 plot, absc, arg, $ xstyle=1, title='phase du gain complexe d antenne' stop endif ; Rangement des harmoniques dans harm_N ; NS0 avec E2, E1, H1 a H16 0 a 17 ; NS1 ---- ---------------- 18 a 35 ; NS8 ---- ---------------- 144 a 161 ; NS23 ---- ---------------- 414 a 431 ; NS45 ---- ---------------- 432 a 449 ; jusqu'au 7 oct 98 inclus: ; E0 ---- NS0 a NS17 450 a 467 ; -- ---- NS10 a NS23 468 a 481 ; -- ---- NS45 482 ; -- ---- NS7 a NS9 483 a 485 ; a partir du 8 oct 98 inclus: ; NS45 ---- NS0 a NS17 450 a 467 ; -- ---- NS10 a NS23 468 a 481 ; -- ---- NS45 482 ; -- ---- NS7 a NS9 483 a 485 ; E2 ---- E2, E1, H1 a H16 486 a 503 ; E1 ---- E2, E1, H1 a H16 504 a 521 ; E0 ---- E2, E1, H1 a H16 522 a 539 ; NS8 ---- NS0 a NS17 540 a 557 ; NS0 ---- NS10 a NS23 et NS45 558 a 572 ; --- ---- NS7 a NS9 573 a 575 ; AA1 ---- NS0 a NS17 576 a 593 ; AA2 ---- ---------- 594 a 611 ; AA3 ---- ---------- 612 a 629 ; AA4 ---- ---------- 630 a 647 ; XX1 ---- YY1 648 a 665 ; XX2 ---- YY2 666 a 683 ; XX3 ---- YY3 684 a 701 ; XX4 ---- YY4 702 a 719 ; Calcul de harm_N produit par le reseau imparfaitement phas'e. ; Rappel du l'utilisation de correl dans RH_MALC_IM_2D : ; bew = xant(correl(1, *)) - xant(correl(0, *)) ; bns = yant(correl(1, *)) - yant(correl(0, *)) ; La base est orientee de l'antenne de reference vers l'antenne de ; reseau. L'indice 0 correspond a la reference, l'indice 1 au reseau. ; La phase phi d'un harm est phase(ant de res) - phase(ant ref), donc: harm_N = harm_N_theo * conj( gain_c(correl(0, *)) ) * $ gain_c(correl(1, *)) ; Controle du facteur de gain (a permis de trouver 0,98 au lieu de 0.98). icon = 0 if (icon eq 1) then begin c0 = gain_c(correl(0, *)) ; antenne reference. c1 = gain_c(correl(1, *)) ; antenne reseau. n1 = 0 & n2 = 36 ; domaine d'indice a tracer. window, 0 & wset, 0 plot, abs(c0(n1:n2)), $ xstyle=1, title='abs(facteur ant de reference)' window, 1 & wset, 1 plot, abs(c1(n1:n2)), $ xstyle=1, title='abs(facteur ant de reseau)' stop endif icon = 0 if (icon eq 1) then begin fac_0 = gain_c(correl(0, *)) fac_1 = gain_c(correl(1, *)) print, 'Avant reform:' help, fac_0, fac_1 fac_0 = reform(fac_0) fac_1 = reform(fac_1) print, 'Apres reform:' help, fac_0, fac_1 plot, abs(fac_0) , title='abs(fac_0) ____, abs(fac_0)... ' oplot, abs(fac_1), linestyle=1 window, 1 & wset, 1 plot, abs(fac_1) , title='abs(fac_1) ____' stop endif ; CONTROLES (AVEC CALCUL PRELIMINAIRES NE SERVANT QU'AUX CONTROLES) ; Calcul de l'image avec harm_N ; freq = 164 ; bidon ie0_ew = 1 ; i_stop = 0 ; 1 pour controles. ; print, 'SIMUL_HARM_N : appel de RH_MALC_IM_2D (pour controles seulement)' if (i_stop eq 1) then print, $ 'SIMUL_HARM_N : appel de RH_MALC_IM_2D (pour controles seulement)' i_recentrage = 0 & x_centrage = 0 & y_centrage = 0 ; bidons. np7 = 128 RH_MALC_IM_2D, $ ; entree i_stop, i_sys_lin, $ i_pave_plein, i_recentrage, x_centrage, y_centrage, $ freq, ie0_ew, harm_N, correl, np7, n_bord, $ ; sortie: image et lobe theorique, centre du Soleil en (np/2, np/2). tf_l_th, l_th, tf_im_2d, im_2d, $ flux_total, flux_compact, uu, vv ; uu, vv coord des frq spatiales dans le plan (u, v). ; print, 'SIMUL_HARM_N : sortie de RH_MALC_IM_2D (pour controles seulement' if (i_stop eq 1) then print, 'SIMUL_HARM_N : sortie de RH_MALC_IM_2D' + $ ' (controles seuls' objet_512 = congrid(objet, 512, 512, cubic=-0.5) im_2d_512 = congrid(im_2d, 512, 512, cubic=-0.5) l_th_512 = congrid(l_th , 512, 512, cubic=-0.5) tf_512 = congrid(tf_im_2d, 512, 512, cubic=-0.5) a_tf_512 = abs(tf_512) ; print, 'SIMUL_HARM_N (sortie) : Flux reconstitues par RH_MALC_IM_2D :' ; print, ' flux total = ', flux_total, $ ; ' (par somme des bordures nulle)' ; print, ' flux compact = ', flux_compact, $ ; ' (par extrapol quadratique de visibilite)' ; Visualisation de l'objet simule et de son image icon = 0 if (icon eq 1) then begin window, 0 & wset, 0 shade_surf, objet_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'SIMUL_HARM_N: objet simule' window, 1 tvscl, objet_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'SIMUL_HARM_N: objet simule' window, 2 shade_surf, im_2d_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'SIMUL_HARM_N: image RH d objet simule' window, 3 tvscl, im_2d_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'SIMUL_HARM_N: image RH d objet simule' stop endif ; Visualisation de l'image de l'objet simule et de sa TF (avec coupe EW) icon = 0 if (icon eq 1) then begin window, 0 & wset, 0 shade_surf, l_th_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'SIMUL_HARM_N: lobe theorique' window, 1 & wset, 1 shade_surf, im_2d_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'SIMUL_HARM_N: image' window, 2 & wset, 2 npz = 52 * 512/np ; nbre de pts zoomes autour de l'origine. shade_surf, a_tf_512(256-npz:256+npz, 256-npz:256+npz) xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'SIMUL_HARM_N: abs(TF(image simulee))' stop window, 0 & wset, 0 h1 = 10 & h2 = 35 absc = h1 + indgen(h2 - h1) plot, absc, abs(tf_im_2d(np/2+h1:np/2+h2, np/2)), xstyle=1 xyouts, /normal, 0.5, 0.97, alignment=0.5, $ ' Coupe EW abs(TF source simulee)' stop endif icon = 0 ; controle des harm de E0 if (icon eq 1) then begin window, 0 & wset, 0 print, 'simul_harm_n : harm_N(523:527) (E0 avec E1 H1 H2 H3 H4)' print, harm_N(523:527) plot, title='(fin simul_* :harm_N(E0)', yrange=[-2, 2], $ abs (harm_N(523:539)) oplot, float (harm_N(523:539)), psym=2 ; 2 *, 3 ., 4 losange oplot, imaginary(harm_N(523:539)), psym=4 window, 1 & wset, 1 tvscl, congrid(abs(visib_rh), 512, 512, cubic=-0.5) stop endif ; Visualisation de l'objet simule et de son image icon = 0 if (icon eq 1) then begin objet_512 = congrid(objet, 512, 512, cubic=-0.5) im_2d_512 = congrid(im_2d, 512, 512, cubic=-0.5) window, 0 & wset, 0 shade_surf, objet_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'simul_harm_N: objet simule' window, 1 tvscl, objet_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'simul_harm_N: objet simule' window, 2 shade_surf, im_2d_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'simul_harm_N: image RH d''objet simule' window, 3 tvscl, im_2d_512 xyouts, /normal, 0.5, 0.95, alignment=0.5, $ 'simul_harm_N: image RH d objet simule' stop endif ; print, 'SIMUL_HARM_N : sortie' end ; fin de SIMUL_HARM_N ;______________________________________________________________________________ PRO SIMUL_RH, $ ; Sortie objet, visib, $ ; objet et sa TF centr'es. l_th, tf_l_th_c, $ ; centres tous les deux. image, tf_im_c ; --------------------- ; But: produire le lobe theorique et l'image d'un soleil a partir d'une repar- ; tition d'harmoniques comportant les antennes anti-aliasing. ; L'image est obtenu par convolution avec un lobe theorique (sans ; erreurs de gain complexe) ; On ne passe plus par les tableaux de positions d'antennes ni par le ; tableau des correlations. print, 'Debut de SIMUL_RH pour simuler : print, ' . un objet soleil et sa TF (avec SIMUL_VISIB_2), print, ' . lobe theorique, image et leurs TF avec les antennes' + $ ' anti-aliasing' ; DEFINITION DU NBRE DE POINTS SUR IMAGE ET SIMULATION DE L'OBJET SIMUL_VISIB_2, $ ; definitions source dans SIMUL_VISIB ; sortie np, $ ; dimension image. objet, $ ; objet simul'e en position relle. visib ; TF (centree dans espace TF) d'objet simul'e ; en position reelle. ; CONSTRUCTION DE LA TF CENTREE DU LOBE THEORIQUE ; Construction de la repartition centree des harmoniques theoriques du RH, ; antennes anti-aliasing exclues tf_l_th_c = complexarr(np, np) pave = complexarr(33, 47) for i = 0, 32, 2 do begin for j = 0, 46 do begin pave(i, j) = 1 endfor endfor tf_l_th_c(np/2-16, np/2-23) = pave ; Ajout des harm 1D EW de E0, E1, E2 tf_l_th_c(*, np/2) = 1 ; Ajout des harm NS de NS45 for j = 0, 45 do begin tf_l_th_c(np/2, np/2-j) = 1 tf_l_th_c(np/2, np/2+j) = 1 endfor ; Ajout des barres EW de NS45 for i = 0, 16 do begin tf_l_th_c(np/2-i, np/2-45) = 1 tf_l_th_c(np/2+i, np/2-45) = 1 tf_l_th_c(np/2-i, np/2+45) = 1 tf_l_th_c(np/2+i, np/2+45) = 1 endfor ; Ajout des barres NS de E2 for j = 0, 23 do begin tf_l_th_c(np/2-48, np/2+j) = 1 tf_l_th_c(np/2+48, np/2-j) = 1 endfor ; Ajout du sous=pave antialiasing avec i0=4 antennes dans le sous-reseau EW ; situe au niveau de l'antenne NS de rang j0 (8 ) ; => harm 1, 3, 5, 7 en EW, ce qui remplit le sous-pave de dimensions ; (2*i0-1)*2 + 1 en EW et 2*j0+1 ; et j0 en NS i0 = 4 j0 = 8 ; rang NS du sous-reseau EW s_pave = complexarr(2*(2*i0-1) + 1 , 2*j0 + 1) s_pave(*, *) = 1 tf_l_th_c(np/2 - (2*i0 - 1), np/2 - j0) = s_pave icon = 0 if(icon eq 1) then begin window, 0 shade_surf, abs(tf_l_th_c) window, 1 tvscl, congrid(abs(tf_l_th_c), 512, 512, cubic=-0.5) stop endif ; FIN DE CONSTRUCTION DE LA TF CENTREE ; TF du lobe theorique centree au coin inferieur gauche tf_l_th = shift(tf_l_th_c, -np/2, -np/2) ; Lobe theorique centree l_th = float(shift(fft(tf_l_th, -1), np/2, np/2)) icon = 0 if(icon eq 1) then begin window, 0 shade_surf, l_th xyouts, 0.5, 0.97, /normal, alignment=0.5, 'Lobe theorique' window, 1 tvscl, congrid(l_th, 512, 512, cubic=-0.5) stop endif ; Calcul de l'image du soleil ; Rappel : les harmoniques acquis par le RH sont echantillonnes sur la ; TF "quasi-reelle" de l'objet. Avec une fft directe (drapeau -1) ; on obtient une image centree au coin inferieur gauche. tf_im_c = visib * tf_l_th_c tf_im = shift(tf_im_c, -np/2, -np/2) image = float(shift(fft(tf_im, -1), np/2, np/2)) ; Visualisation de l'objet simule et de son image icon = 0 if(icon eq 1) then begin objet_512 = congrid(objet, 512, 512, cubic=-0.5) image_512 = congrid(image, 512, 512, cubic=-0.5) window, 0 shade_surf, objet_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'Objet simule' window, 1 tvscl, objet_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'Objet simule' window, 2 shade_surf, image_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image avec sous-pave anti-aliasing' window, 3 tvscl, image_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image avec sous-pave anti-aliasing' stop endif end ; fin de SIMUL_RH ;______________________________________________________________________________ PRO RH_SIMULATION ; mettre en commentaire pour demarrer directement ; (et appeler ensuite simul_rh) ; But : procedure bidon, appelee dans le seul but de compiler les procedures ; contenues dans le fichier. J_M_Delouis = 0 ; in memoriam. end ; fin de RH_SIMULATION. ;______________________________________________________________________________