PRO FILTRE_PAVE, image, je0_ew, $ ; entree image_2 ; sortie ; But : retirer de l'image les contributions des harmoniques depassant du pave ; et les harm de E0 interieurs au pave si je0_ew=1 . ; Ceci dans l'espoir insense de calculer le flux en supposant que la ; somme S et N des bordures ded image_2 est nulle. taille = size(image) np = taille(1) tf = fft(image, 1) ; origine de TF au coin inf gauche. tf_c = shift(tf, np/2, np/2) ; Calcul du filtre centre filtre = fltarr(np, np) pave = replicate(1, 33, 47) if (je0_ew eq 0) then begin for i=1, 31, 2 do begin pave(i, 23) = 0 ; for i=0, 15 do begin ; pave(2*i+1, 23) = 0 endfor endif filtre(np/2-16, np/2-23) = pave icon = 0 if (icon eq 1) then begin f_512 = congrid(filtre, 512, 512, cubic=-0.5) window, 0 tvscl, f_512 window, 1 shade_surf, f_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'filtre du pave' stop endif ; Filtrage de l'image tf_c_f = tf_c * filtre tf_f = shift(tf_c_f, -np/2, -np/2) image_2 = float(fft(tf_f, -1)) end ; fin de FILTRE_PAVE. ;______________________________________________________________________________ PRO RH_MALC_IM_2D, $ ; entree : i_stop, i_sys_lin, $ ch_polar, $ i_pave_plein, i_recentrage, k_max, l_max, $ freq, ie0_ew, harm_N, correl, npi, n_bord, $ ; sortie : k_max_auto, l_max_auto, $ harl, l_th, harm, im_int, $ flux_total, flux_compact, uu, vv ; Procedures appelantes : ; - RH_CALC_IMAGE_HELIO (une fois) ; - SIMUL_HARM_N (dans ~/calib/rh_simulation) pour des controles ; - CALIBRATION (une fois) ; - CALIB_CON ------------------ ; - CALIB_SUB ; - GET_DATA_UV_RH (dans gmrt) ; But : calcul de l'image interferometrique. ; ETAPES DE RH_MALC_IM_2D ; - remplissage du plan uv ; - 1ere reconstitution des composantes continues d'image et lobe theorique ; par extrapolation. ; - calcul d'image interferometrique (avec flux = total des points). ; - amelioration de la composante continue (total des bordures N et S nulle). ; Modifications: ; 99 mai 6: Reconstitution approximative de la composante continue, pour l'ima- ; ge et le lobe theorique, comme moyenne des bas harmoniques. Cela ; permet, dans dpatchfits, de normaliser l'image a celle obtenue ; avec un lobe d'integrale 1 en utilisant directement le lobe theo- ; rique (auparavant depourvu de composante continue) au lieu de son ; fit par un gaussienne, utilise dans la version J-M parce qu'il a ; necessairement une CC. ; 27: Mise en coherence de tout le calcul d'image: ; - definition de positions des antennes dans un repere (inverse ; pour le terrain du dessus) d'axes vers W et N (et non plus W ; et S), identique a celui utilise pour l'image du soleil. ; - definition du vecteur de la base d'un couple d'antennes oriente ; de l'antenne de reference (NS) vers l'autre antenne (les defi- ; nitions etaient opposees pour bew et bns dans la version J-M). ; - remplissage du tableau des harm a 2 dim: la visibilite du couple ; de base (bew, bns) est attribue a l'harmonique (bew, bns) et ; non plus au couple (-bew, bns) comme dans la version J-M). ; - usage d'une FFT directe au sens d'IDL (drapeau -1) pour calculer ; l'image (et non plus d'une FFT inverse de drapeau +1 dans la ; version J-M). Ceci parce que bew et bns interviennent avec le ; signe + dans le calcul theorique de phi (ant de ref omega*t, ; 2eme ant omega*t+phi). ; 99 sep 7: normalisation de facon que la somme des pts de l'image interfero- ; metrique soit le flux. ; 00 avr 10: Determination du flux en imposant que la somme des points de ; l'image aux lisieres nord et sud (jamais affectees par ; l'aliasing) soit nulle. ; nov 30: Deplacement de l'appel a REGCLEAN a l'exterieur de MALC_IM_2D, ; dans CALC_IMAGE_HELIO, juste apres l'appel a MALC_IM_2DS. ; Supression des arguments et mots-cle de la liste de MALC_IM_2D ; qui ne servaient qu'a REGCLEAN. ; 01 mar 19: Ajustement du flux avec la somme des bordures N et S nulles. Les ; harm de E0 donnent une contribution qui ne s'annule pas sur ; ces bordures, => si E0 est dephase le flux peut etre tres faux. ; mai 31: Introduction de i_sys_lin (passe par CALIBRATION et transmis a ; FIT_QUADRIC), qui permet de choisir la methode de resolution ; des systemes lineaires (Cramer, LU, SVD). ; 01 jul 10: npew est change en npi (dimension de l'image interferometrique) ; qui devient distinct de la dimension np de l'image heliographi- ; que. Mise a jour des notations dans MALC_IM_2D et HARM_N_VIS. ; 17: Commentaires sur l'utilisation de fft directe (drapeau -1) et ; inverse (drapeau 1) pour les images du RH. Correction de qques ; commentaires faux. ; 01 sep 12: Correction : on calculait la base en annulant les bordures N et S ; de l'image mais on n'utilisait pas le resultat ! ; La hauteur n_bord de ces bandes est passee en argument ; 13: La definition de n_bord, pour ajuster la base dans MALC_IN_2D, ; est remontee dans DPATCHFITS et transmise aux procedures appe- ; lantes (CALIBTATION, CALC_IMAGE_2D, SIMUL_HARM_B, ; CALIB_CON_10) ; 20: 2 arguments supplementaires en sortie : ; . le flux total, composante continue obtenue par la methode de ; somme nulle sur les bordures N et S ; . le flux compact, composante continue obtenue par extrapola- ; tion quadratique de la visibilite a l'origine. ; 25: complement du pave central par interpolation (pour reduire l'alia- ; sing). ; 01 oct 3: restauration des valeurs initiales de harm sur ,l'axe u apres le ; completement du pave (ces harm impairs sont connus et n'ont ; pas besoin d'etre interpoles. ; 11: possibilte de substituer une distribution d'harmoniques correspon- ; dant a des antennes anti-alias (et un reseau parfaitemment ; phas'e) a la distribution obtenue avec harm_N. ; 03 mar 12: amelioration du remplissage du pave central : on interpole la TF ; de l'image dont le max a ete place au centre des coordonnees. ; On calcul ensuite la TF de l'image remise en place, avec le ; centre du soleil au centre des coordonnees. ; ajout de harm et harl a la liste de sortie. ; ajout de i_recentrage a la liste d'entree. ; 03 mar 19: ajout de k_max et l_max (appeles x_centrage et y_centrage dans ; les procedures en amont) a la liste d'entree, ; i_recentrage peut etre egal a 0 (pas de recentrage), 1 (recentrage ; sur le max d'image), 2 (recentrage sur les coordonnees inter- ; ferometriques x_centrage et y_centrage (par rapport au centre ; du champ). ; 03 dec 26: extension aux antennes AA, en gardant la compatibilite avec les ; observations d'avant les antennes AA. ; 05 dec 2: amelioration de l'estimation de la base de l'image : la base est ; d'abord estimee en annulant la somme des points de l'image dans ; 2 bandes N et S, de largeur (choisie dans RH_DPATCHFITS_NRH) de ; 10 canaux interferometriques (sur 128). Mais une evaluation ; imprecise de "moy" peut fortement affecter le flux, obtenu par ; integration sur tout le champ interferomerique CI, puisque ; l'erreur sur le flux resulte de l'integration sur tout le CI ; CI, alors que le soleil est nettement moins etendu (c'est enco- ; re plus vrai pour un centre d'orage tres compact). ; Pour reduire cet effet, ayant deja obtenu une 1ere approximation ; du niveau zero en prenant nulle la moyenne de l'image sur les ; bandes N et S, on integre le flux de la source sur un domaine ; ou la source est notable (p ex. > 10% de son maximum par rap- ; port au niveau zero approximativement determin'e), mais sans ; exclure les secondaires negatifs proches du max d'une source ; compacte, ni les secondaires >0 ou <0 un peu lointains mais ; non negligeables (qques %). Comentaires plus bas (rechercher ; la chaine "sur un domaine ou l'image est notable" ; Ainsi "flux_total" ne resulte plus d'une addition sur tous les ; points de l'image interferometrique et n'est pas (fortement) ; affecte par une petite imprecision sur le niveau zero obtenu ; par la methode des bandes N et S. ; 06 mar 2: le recentrage pour remplir le pave central est fait sur sur le ; maximum de la valeur absolue d'image, et non plus sur le maxi- ; mum d'image, pour pouvoir traiter une polar negative. ; chgt de "image" en "im_int" (comme image interferometrique). ; 06 mar 8: le recentrage automatique polar est pris identique au recentrage ; automatique polar, considere comme plus fiable (le 14 nov 2006 ; le centrage auto polar etait faux de 1/2 periode en EW => echan ; ge de l'alias et de la source) ; ajout de ch_polar, x_centr_auto, y_centr_auto dans la liste. ; Notations: ; ch_polar 'non polar' ou 'polar'. ; i_pave_plein 1 pour remplir le pave central par interpolation. ; i_recentrage recentrage d'image avant interpolation de TF, ; 0 sans recentrage ; 1 avec ---------- sur son maximum de valeur absolue d'image ; 2 ----------------sur coord interfero ci-dessous. ; x_centrage abscisse de recentrage choisie en Rs dans RH_DPNEW et trans- ; formee en en canaux interfero (de -npi/2 a + npi/2 - 1) par ; RH_DPATCHFITS_NRH. Sert dans la procedure de remplissage du ; pave central. ; Pour le recentrage auto polar, x_centrage transmet la posi- ; tion de recentrage automatique calculee pour le non polar, ; considere comme plus fiable. ; k_max absc interferometrique en canaux p r centre de l'image. ; l_max ord -------------------------------------------- ; nbre_harm nombre d'harmoniques dans une scrutation (576 ou 648 ou 720). ; harm_N complexarr (nbre_harm) harmoniques, rang'es au format Nancay. ; Rempli en non polar ou polar juste avant l'appel a ; RH_CALC_IMAGE_HELIO. ; harm, harl complexarr (128, 128) tableaux des harm dans le plan (u, v). ; En sortie ils sont centres en milieu de tableau. ; correl tableau (2, nbre_harm) indiquant les numeros des 2 antennes ; correspondant a chacun des nbre_harm harmoniques. Permet de ; faire le passage de harm_N a harm. ; freq frequence (MHz). Ne sert que pour mise au point (flux du Cygne) ; im_int image interferometrique brute de npi*npi points, normalisee de ; sorte que son flux soit le total de ses points). ; i_stop 1 pour controles et arrets. ; l_th lobe theorique, normalise a un total 1 (avec un max de l'ordre ; de 0.06). ; C'est l'image (avec composante continue) d'une source dont ; l'amplitude des harmoniques est 1, c'est a dire, compte tenu ; de la procedure de calibration, d'une source ponctuelle de ; 1 sfu. ; Rem : cette normalisation est differente de celle du lobe theo- ; rique filtre d'echelle l_th_ech qui sera normalise a un ; maximum de 1 dans WCLEAN, ce qui correspond a l'usage ; qui en est fait dans P_CLEAN_2D pour soustraire des com- ; posantes clean de l'image residuelle. ; npi dimension dee l'image interferometrique brute (non cleanee) ; "im_int". Depuis le 6 jul 01 npi est fige a 128 et distinct ; de la dimension np de l'image heliographique. ; x_centr_auto position de centrage automatique calculee dans le cas non polar ; y_centr_auto pour etre fournis au cas polar sous les noms de x_centrage ; et y_centrage. ; uu, vv coord de frequences spatiales dans le plan (u, v). ; print, 'RH_MALC_IM_2D entree : k_max, l_max =', k_max, l_max ;_____________________________________________________________________________ ; HARMONIQUES DU RH, FFT EN IDL, CALCUL DE L'IMAGE : POSITION DU PROBLEME ; 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 a un signe - dans ; l'exponentielle complexe. ; - ----- inverse (------ +1) ne comporte pas ------- 1/n et a un ---- + dans ; ------------------------. ; - Proprietes : ; 1) TF inverse (0) = total( fonction de depart) ; 2) fonction de depart (0) = total (TF directe) ; 3) TF directe(TF invserse) = TF invserse (TF directe) = identite ; 4) TF directe(TF directe) = 1/n * identite ; 5) TF invserse(TF invserse) = n * identite ; Rappel sur les harm obtenus par le RH (reprises en partie de RH_HARM_N_VIS) : ; a) correl(0, i) et correl(1, i) sont les numeros des antennes de reference ; et de reseau. ; b) les bases bew et bns, en unites des pas des reseaux dew et dns, sont ori- ; entees de l'antenne de reference vers l'antenne de reseau (modification ; du 27 mai 99). ; c) l'amplitude de la composante continue et des bas harmoniques est egale ; au flux (sfu). ; ; Precisions sur le formalisme. On appelle : ; . Wew et Wns les vecteurs unitaires portes par les reseaux, approximative- ; ment vers W et N. ; . phi_a le dephasage aerien pour un coupe d'antenne (omega*t a l'antenne ; de reference, omega*t + phi_a a l'antenne de reseau) ; On a : ; phi_a = 2pi/lambda (bew dew Wew.u + bns dns Wns.u) (1) ; ou dew et dns sont les pas EW et NS des reseaux, ; bew et bns -------- bases normalisees (entieres) ; On suppose, comme on l'a toujours fait, que la sortie du correlateur comple- ; xe fait intervenir exp(i*phi_a) . ; Explicitons les frequences spatiales. On pose ; u = u0 + u1 (u0 vers le centre de la source) ; et on somme sur la source etendue. La sortie S du correlateur est : ; S = exp(i 2pi/lambda D.u0) Somme[B(u1) * exp(i 2pi / lambda * ; (bew dew Wew.u1 + bns dns Wns.u1) domega)] (2) ; ou B est la brillance par rapport a 2 coordonnes angulaires, ; domega l'angle solide differentiel, ; bew et bns les bases EW et NS non normalisses par les pas des reseaux. ; On definit les "coordonnees interferometriques" X et Y ; X = dew Wew . u1 / lambda et Y = dns Xns . u1 / lambda ; et une nouvelle brillance Bi par rapport a X et Y ; B(u1) domega = Bi(X, Y)dX dY ; L'equation (2) s'ecrit, en faisant disparaitre la poursuite des franges : ; S = Somme[Bi(X, Y) * exp(i 2pi / lambda (bew X + bns Y)) dX dY] (3) ; Les frequences spatiales conjugees de X et Y, multiples de l'unite, sont : ; uu = bew et vv = bns (4) ; ; Indication de P. Picard le 30 avril 03 : dans un couple d'antennes, l'an- ; tenne de reference est celle qui porte le dipole destine a la polar. ; Elle est reliee a l'entree du correlateur qui comporte un retard sup- ; plementaire 1/4 de periode. Cela fait, a 1ere vue, que la sortie du ; correlateur complexe comprte exp(-i*phi_a), et est donc le complexe ; conjugue de ce qu'on a toujours suppose. ; Possibilites : ; . une erreur sur la correction oscillateur inferieur/superieur, qui ; remettrait les choses a l'endroit, ; . c'est peut-etre pour ca qu'il a fallu remplir les tables de gain ; avec conj(gains) et non les gains. ; Conclusions pour le calcul de l'image : ; Si on laisse tomber la remarque de P. Picard et qu'on utilise : ; - la propriete 1) des fft IDL ; TF inverse (0) = total( fonction de depart) ; jointe a la propriete c) xdes harmoniques fournis par le RH ; - l'expression (3) de la sortie du correlateur complexe ; on conclut que ; - la TF inverse (drapeau 1) realise la meme operation que celle faite par ; le recepteur du RH. ; - la TF directe (drapeau -1) calcule l'image "interferometrique" de la sour- ; ce, normalisee de sorte que la somme des points soit egale a son flux ; (sfu), a partir des harmoniques fournis par le RH (propriete 3). ; En pratique on ne calcule le flux que quand l'image est exprimee en en sfu ; (et non en Kelvins). ; ; RH_HARM_N_VIS remplit le plan uv selon les relations (4). ;______________________________________________________________________________ ; DEBUT DE LA PARTIE ACTIVE de RH_MALC_IM_2D version_solarsoft = 1 ; 0 version test de Claude 1 SolarSoft. ; permet de gerer impressions et stop nbre_harm = n_elements(harm_N) ; Controle des harm EW de E0 (recopiee de simul_harm_N) icon = 0 ; controle des harm de E0. Trouve faute de signe dans ; RH_SIMUL_HARM_N le 10 fev 2000. if ((icon eq 1) and (i_stop eq 1)) then begin window, 0 & wset, 0 plot, title='RH_MALC_IM_2D: 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=3 print, 'MALC_IM_2D entree: harm_N(523:527) (E0 avec E1 H1 H2 H3 H4' print, harm_N(523:527) stop endif ; print, 'RH_MALC_IM_2D : entree :' ; Remplissage du plan uv avec les harmoniques observes RH_HARM_N_VIS, $ ; fichier RH_HARM_N_VIS. i_stop, $ ; entree npi, correl, harm_N, $ ; entree harm, harl, uu, vv ; sortie ; Produit les repartitions "harm" et "harl" des harmoniques spatiaux ; "interferometriques" 2D (repartis sur une grille a mailles carres) ; corrigees des redondances et CENTRES EN MILIEU DE TABLEAU, respec- ; tivement pour l'image (harm) et le lobe theorique (harl). ; Processus de remplissage : on cumule, en chaque point du plan uv (de ; coordonnees u, v), la valeur de l'harmonique u=bew, v=bns. Simul- ; tanement on cumule la valeur conjuguee en u=-bew, v=-bns. La cou- ; verture complete utilisant la symetrie hermittique se fait donc ; harmonique par harmonique, et non pas en symetrisant le 1/2 plan ; directement fourni par les lignes de bases du "T". ; Cela permet d'eviter une difficulte pour les points du pave sur ; l'axe u, pour lesquels il y a redondance des harmoniques de E0 ; et de NS0_ns (ces derniers eux-memes redondants). ; ; Remarque sur les harm EW purs de NS0_ns avec les moities E et W du ; reseau. On appelle : ; . phi_a 2pi*bew/lambda (c'est la phase aerienne d'une antenne ; de reseau de la moitie ouest du reseau EW). ; . phi_r la phase d'une antenne H (de signe oppose a son exces ; de longueur de cable. ; . phi_0 ---------de NS0_ns -------------------------- ; -------------------- ; Les sorties de correlateur d'un harmonique de NS0_ns avec une H de ; la moitie W et la H symetrique de la moitie E sont respective- ; ment : ; Sw = exp i*[( phi_a + phi_r_w - phi_0)] ; Se = exp i*[(-phi_a + phi_r_e - phi_0)] ; Le resultat cumul'e du remplissage du plan uv au point correspon- ; dant a la H de la moitie W (sur le 1/2 axe u>0) est donc : ; Vis = Sw + conj(Se) = exp[i*phi_a] * ; [ exp(i( phi_r_w - phi_0)) ; +exp(i(-phi_r_e + phi_0)] ; Le crochet se met sous une forme plus claire en posant : ; p1 = phi_r_w - phi_0 ; p2 = -phi_r_e + phi_0 ; puis : ; p_S = (p1 + p2) / 2 comme "somme" ; p_D = (p1 - p2) / 2 ----- "difference" ; On a alors : ; Vis = exp[i*phi_a] * [exp(i(p_S + p_D)) + exp(i(p_S - p_D)) ; = -------------- exp(i*p_S) * [exp(i*p_D) + exp(-i*p_D)] ; ; = 2 * exp[i*phi_a] * exp[i*(phi_r_w - phi_r_e)] * ; cos[phi_0 - (phi_r_w + phi_r_e)/2] ; On remarque que : ; . la phase de NS0_ns disparait du facteur de phase (ca tombe ; bien pour NS0_ns puisqu'on 'a choisie comme origine des ; phases du fait de sa position centrale, ca tombera moins ; bien pour NS8) mais subsiste dans le facteur reel d'ampli- ; tude, ; . les phases des antennes de reseau interviennent par leur dif- ; ference dans le facteur de phase. ; ; La redondance de ces harmoniques est corrigee. ; Controle de redondance. Recherche du creux suivant axe v ; c'etait du a une faute de la correction des harm de E0 dans ; CORR_GAIN_ANT, qui affectait les 8 premiers harmoniques NS de NS8. icon = 0 if ( (i_stop eq 1) and (icon eq 1) ) then begin taille = size(harm) n_dim2 = taille(2) / 2 n_z = 32 n1z = n_dim2 - n_z & n2z = n_dim2 + n_z - 1 ;; toto = harm (n1z : n2z, n1z : n2z) toto_512 = congrid(abs(toto), 512, 512) window, 0 loadct, 13 tvscl, toto_512 xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'abs(visibilite acquise) dans RH_MALC_IM_2D' stop ; Le 21 fev 2006, on constate que la redondance est bien corrigee ; et qu'il n'y a pas ce creux selon l'axe v. endif ; Controle des harm de bases 50m icon = 0 ; supression des autres verif apres la faute trouvee le ; 10 fev 2000. if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_MALC_IM_2D verif 1: TF avant calcul d image' nh = 8 ; zoom de -nh a +nh. absc = indgen(2*nh + 1) - nh h_ew_red = harm(npi/2 - nh : npi/2 + nh, npi/2) h_ns_red = reform(harm(npi/2, npi/2 - nh : npi/2 + nh)) a_h_ew_red = abs(h_ew_red) a_h_ns_red = abs(h_ns_red) p_h_ew_red = 180/!pi * $ atan(imaginary(h_ew_red), float(h_ew_red + 0.0001)) p_h_ns_red = 180/!pi * $ atan(imaginary(h_ns_red), float(h_ns_red + 0.0001)) ; Trace des modules des coupes EW et NS window, 0 amax = max(a_h_ew_red) > max(a_h_ns_red) amin = min(a_h_ew_red) < min(a_h_ns_red) da = amax - amin ymin = amin - 0.1 * da & ymax = amax + 0.1 * da plot, absc, a_h_ew_red, xstyle=1, yrange=[ymin, ymax], ystyle=1 oplot, absc, a_h_ns_red, linestyle=1 ; psym = 2 *, 3 ., 4 losange xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'RH_MALC_IM_2D: abs(harm) coupe EW ___, coupe NS ...' ; Trace des phases des coupes EW et NS window, 1 pmax = max(p_h_ew_red) > max(p_h_ns_red) pmin = min(p_h_ew_red) < min(p_h_ns_red) dp = pmax - pmin ymin = pmin - 0.1 * dp & ymax = pmax + 0.1 * dp plot, absc, p_h_ew_red, xstyle=1, yrange=[ymin, ymax], ystyle=1 oplot, absc, p_h_ns_red, linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'RH_MALC_IM_2D: phase(harm) coupe EW ___, coupe NS ..., deg.' stop endif ; fin de controle ; Remplacement eventuel de la distribution d'harmoniques dans le plan (u,v) ; par une distribution avec antennes anti-aliasing ; Cela ne peut se faire qu'ici car la distribution actuelle (sans les ; antennes anti_aliasing) est gravee au burin dans la structure de ; harm_N et les valeurs de "correl". ; RESERV'E A LA PRODUCTION DE TRANSPARENTS SIMULANT DES IMAGES AVEC ; LES ANTENNES ANTI-ALIAS. i_sim_anti_alias = 0 if (i_sim_anti_alias eq 1) then begin ; ecrasement avec SIMUL_RH ; print, 'Appel de SIMUL_RH pour simuler : ; print, ' . un objet soleil et sa TF (avec SIMUL_VISIB_2), ; print, ' . lobe theorique, image et leurs TF avec lea antennes' + $ ; ' anti-aliasing' RH_SIMULATION ; au cas ou on ne l'aurait pas appele avant. SIMUL_RH, $ objet_simule, visib, $ ; objet et sa TF centree. l_th_0 , harl, $ ; centres tous les deux. im_2d_0, harm ; --------------------- ; Rem : l_th_0 et im_2d_0 sont fournis de toutes facons en ; sortie. Il serviront peut-etre a des verifications. flux_total_0 = total(im_2d_0) ; pour verification. taille = size(im_2d_0) np = taille(1) if (np ne npi) then begin print, 'MALC_IM_2D : npi et np differents' stop endif harm(np/2, np/2) = 0 ; comme en sortie de RH_HARM_N_VIS. harl(np/2, np/2) = 0 ; -------------------------------- endif ; Controle amplitude de harm et harl icon = 0 if ((icon eq 1) and (i_stop eq 1) and (i_sim_anti_alias eq 1)) then begin print, 'RH_MALC_IM_2D : i_sim_anti_alias =', i_sim_anti_alias nz = 20 ; zooming de -nz a + nz. h_l = harl(np/2-nz:np/2+nz, np/2-nz:np/2+nz) h_m = harm(np/2-nz:np/2+nz, np/2-nz:np/2+nz) harl_512 = abs(congrid(h_l, 512, 512, cubic=-0.5)) harm_512 = abs(congrid(h_m, 512, 512, cubic=-0.5)) window, 0 tvscl, harl_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(TF(l_th) avec sous-pave anti-aliasing' window, 1 tvscl, harm_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(TF(mage)) avec sous-pave anti_aliasing' stop endif ; Reconstitution du flux compact (composante continue d'image et du lobe ; theorique) par : ; - extrapolation d'ordre 0 (moy des harm voisins de O). 6 mai 99. ; - ------------- parabolique anisotrope des bas harm. 20 sep 99. ; Doit se faire apres la compensation de redondance (dans harm_N_vis). ; Une meilleure determination de la composante continue, incluant le ; flux des grandes structures est faite plus bas an assujetissant la ; moyenne des points aux bords N et S de l'image (jamais affectes par ; l'aliasing) a etre nulle. harl(npi/2, npi/2) = 1 ; simple pour le lobe theorique ! ; Extrapolation d'ordre 0 (en EW les harm impairs sont nuls hors axe sans les ; antennes AA. Je ne retouche pas le programme pour les antennes AA : ; l'extrapolation est simplement faite sur un peu moins de points qu'on ; pourrait le faire). if (ie0_ew eq 0) then begin ; 8 harm: -2 0 2 en EW, -1 0 1 en NS. c_c_0 = $ ( abs(harm(npi/2+2, npi/2 )) + abs(harm(npi/2+2, npi/2+1)) + $ abs(harm(npi/2 , npi/2+1)) + abs(harm(npi/2-2, npi/2+1)) + $ abs(harm(npi/2-2, npi/2 )) + abs(harm(npi/2-2, npi/2-1)) + $ abs(harm(npi/2 , npi/2-1)) + abs(harm(npi/2+2, npi/2-1)) ) / 8 endif if (ie0_ew eq 1) then begin ; 4 harm: -1 et 1 en EW, -1 et 1 en NS. c_c_0 = $ ( abs(harm(npi/2+1, npi/2 )) + abs(harm(npi/2 , npi/2+1)) + $ abs(harm(npi/2-1, npi/2 )) + abs(harm(npi/2 , npi/2-1)) ) / 4 endif harm(npi/2, npi/2) = c_c_0 ; eventuellementb ecrase par c_c_p. flux_compact = c_c_0 ; -------------------------------- ; Extrapolation parabolique anisotrope. i_extrapol_quad = 0 ; 1 (en principe meilleure evaluation du flux ; compact peut quelqufois planter. if (i_extrapol_quad eq 1) then begin n_harm = 10 ; nbre d'harm de part et d'autre de l'origine sur lesquels on fait ; le fit parabolique de la visibilite. ; Le fit est reduit dans RH_FIT_QUADRIC aux pts > fraction (0.2) ; de la moyenne de points autour de l'origine, voir details dans ; RH_FIT_QUADRIC. Cette fraction est ramenee de 0.4 a 0.2 le 23 ; oct 02 dans l'espoir que ca plante moins souvent dans ; RH_LINEAR_LU. Sinon sauter carrement l'extrapolation quadrati- ; que. ; Rem : Ca peut planter dans CRAMER appele par RH_FIT_QUADRIC ; (input array is singular). Reduire n_harm. j_stop = i_stop ; i_stop deja utilise en entree ; pour eviter arret non desire a l'appel de MALC_IM_2D pour ; controles dans SIMUL_HARM_N. RH_FIT_QUADRIC, j_stop, i_sys_lin, $ ; entree abs(harm), n_harm, $ ; entree tcon, tnul, n_tnul, $ ; sortie. con pour "contrainte" c_c_p, BB, CC, DD, fit_visib ; sortie ; Calcul de la cc par un fit parabolique isotrope a l'origine : ; fit_visib = c_c_p + BB*i^2 + CC*i*j + DD*j^2 ; BB, CC et DD ne sont pas utilises ici. ; Le mot cle ident ne figure pas (utilise seulement pour le cas ; 1D). ; Les variables tcon, tnul et n_tnul sont produites en sortie pour ; RH_FIT_CUBIC pour le fit d'une phase de visibilite quand on a ; calcul'e le fit du module de cette visibilite avec RH_FIT_- ; QUADRIC (en fait on utilise plutot RH_FIT_QUARTIC et RH_FIT_- ; QUADRIC n'est plus utilise qu'ici pour calculer la composante ; continue). if (i_stop eq 1) then begin ; essais seulement. ; pour eviter ce qui suit dans l'appel de RH_MALC_IM_2D par ; SIMUL_HARM_N, avant lequel il faut alors faire i_stop=0. ; c_c_p = 0 ; essais seulement ; if (c_c_p eq 0) then print, $ ; 'WARNING : c_c_p mis a zero dans RH_MALC_IM_2D' endif harm(npi/2, npi/2) = c_c_p ; ecrase eventuellement c_c_0 . flux_compact = c_c_p endif ; fin de l'extrapolation parabolique. ; Choix entre les 2 calculs de la composante continue. ; FIN DE RECONSTITUTION DU FLUX COMPACT (EXTRAPOLATION DE LA COMP CONTINUE) ; Controle de la reconstitution de la composante continue icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin harm_c = abs(harm(npi/2 - n_harm : npi/2 + n_harm, $ npi/2 - n_harm : npi/2 +n_harm) ) ; Rappel : n_harm est le nbre d'harm de part et d'autre de ; l'origine sur lequel on extrapole la comp continue. dom_0 = where(harm_c eq 0) ; help, harm_c ; Visualisation de la partie centrale du module de la visibilite et compa- ;raison avec le fit fit_2 = fit_visib & fit_2(dom_0) = 0 ; pour visu de la diff harm_c_512 = congrid(harm_c, 512, 512, cubic=-0.5) fit_512 = congrid(fit_2 , 512, 512, cubic=-0.5) window, 0 shade_surf, harm_c_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'abs(visibilite observee)' window, 1 shade_surf, fit_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'fit' window, 2 shade_surf, fit_512 - harm_c_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'fit - abs(visib observee)' ; Impression du flux ch_format = "('composante continue extrapolation ordre 0 =', " + $ "f8.3, ' , ordre 2 =', f8.3)" print, format=ch_format, c_c_0, c_c_p stop endif ; Fin de controle. ; Controle : visualisation de l'amplitude des harmoniques (histogramme) icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin wset, 0 ; Zoom sur la partie centrale des harmoniques (de -nh a + nh). nh = 18 h_i_c = harm(npi/2-nh : npi/2+nh, npi/2-nh : npi/2+nh) ; pour image h_l_c = harl(npi/2-nh : npi/2+nh, npi/2-nh : npi/2+nh) ; pour lobe ; Trace du tableau des harmoniques et histogramme. amul = (freq / 164)^0.9 ; pour corriger le spectre du Cygne. tvscl, congrid(amul * abs(h_i_c), 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harmoniques observes) (impairs interpoles)' window, 1 & wset, 1 bin_size=0.01 plot, histogram(amul * abs(h_i_c), binsize=bin_size), $ xrange=[0.1/bin_size, 1.5/bin_size] ; stop ; plot, h_i_c(*, 64 + 46) ; stop ; Idem sur le lobe theorique. ; window, 2 & wset, 2 tvscl, congrid(abs(h_l_c), 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(TF lobe theorique) (harm impairs interpoles)' tab_z = where(abs(harm_N) eq 0., c_z) ; c nbre d'harm nuls. if (c_z lt 1) then print, 'nbre d harm nuls: 0' if (c_z ge 1) then begin print, 'nbre d harm nuls:', c_z, ' , de numeros:' print, tab_z endif stop ; Conclusions: ; - Les correlations de E0 et des NS sont bien presentes dans ; l'observation du 25 jul 98. Sont nulles les sorties corres- ; pondant aux correlations: ; . E_NS (NS45) avec E2 a H16 ; . E0 avec NS45 ; . E2 avec H1 (panne?) ; . E1 avec E1 ; . NS8 avec NS8 ; . NS0 avec E_NS ; - Les amplitudes des harmoniques du Cygne sont de l'ordre de ; 164/freq pres de l'origine (Cygne du 16 jul 98). endif ; Fin de mise au point ; DEBUT DE CALCUL ET TRAITEMENT D'IMAGE AVEC PAVE LACUNAIRE, OU PARTIELLEMENT ; LACUNAIRE SI ON A LES ANTENNES AA (apres le 26 dec 03) ; . calcul de l'image, ; . ajustement de la composante continue par la methode des bandes. ; Controle de redondance. Recherche du creux suivant axe v ; c'etait du a une faute de la correction des harm de E0 dans ; CORR_GAIN_ANT, qui affectait les 8 premiers harmoniques NS de NS8. icon = 0 if ( (i_stop eq 1) and (icon eq 1) ) then begin taille = size(harm) n_dim2 = taille(2) / 2 n_z = 32 n1z = n_dim2 - n_z & n2z = n_dim2 + n_z - 1 toto = harm (n1z : n2z, n1z : n2z) toto_512 = congrid(abs(toto), 512, 512) window, 0 loadct, 13 tvscl, toto_512 xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'abs(visibilite acquise) dans RH_MALC_IM_2D' stop ; Le 21 fev 2006, on constate que la redondance est bien corrigee ; et qu'il n'y a pas ce creux selon l'axe v. ; On constate aussi que la composante continue a ete reconstituee. endif ; Calcul d'image et de lobe theorique avec centre du Soleil en (npi/2, npi/2) im_int = float (shift(fft(shift(harm, -npi/2, -npi/2), -1), npi/2, npi/2)) l_th = float (shift(fft(shift(harl, -npi/2, -npi/2), -1), npi/2, npi/2)) ; Rappel : FFT directe (drapeau -1) pour obtenir une image. ; Rem : total(im_int) = comp continue d'apres la propriete 2 des fft en ; IDL. Si la comp continue utilisee est mauvaise l'image est ; simplement decalee. ; Controle image finale et lobe avant ajustement de la composante continue icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin i_sh_tv = 2 ; 1 shade_surf, 2 tvscl ; Trace du lobe if (i_sh_tv eq 1) then begin window, 0 shade_surf, congrid (l_th, 512, 512, cubic=-0.5) endif if (i_sh_tv eq 2) then begin window, 0, xsize=512, ysize=512 tvscl, congrid (l_th, 512, 512, cubic=-0.5) endif xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'Lobe theorique' ; Trace de l'image if (i_sh_tv eq 1) then begin window, 1 shade_surf, congrid (im_int, 512, 512, cubic=-0.5) endif if (i_sh_tv eq 2) then begin window, 1, xsize=512, ysize=512 tvscl, congrid (im_int, 512, 512, cubic=-0.5) endif xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'Image avant ajustement de la composante continue' stop endif ; Fin controle image et lobe avant ajustement de la composante ; continue. ; Calcul du flux total par somme des bordures N et S nulles (image avec pave ; lacunaire) ; Essais du 19 fev 01 a 19:45 (avant filtrage inclus dans malc_im_2d): ; nz = 10 Flux par bordure nulle ; Comp large seule (flux=1) ; phas'e non phas'e ; ie0_ew = 0 flux = -0.56 -0.47 ; ie0_ew = 1 1.007 0.92 ; Comp compacte seule (flux=1) ; phas'e non phas'e ; ie0_ew = 0 flux = -2.22 -2.12 ; ie0_ew = 1 1.016 0.74 ; On ajoute une constante a l'image de facon que la moyenne des points ; sur 2 lignes au bord N et 2 lignes au bord S de l'image (bords qui ; ne jamais affectes par l'aliasing) soit nulle. ; Essais du 23 mars 01 : ; - comp compacte : flux_quad flux_nf flux_f avecE0 sans E0 ; . phas'e 1.004 0.996 1.029 -1.07 ; . E0 seul phas'e 1.025 1.064 1.110 -1.03 ; . non phas'e 1.019 0.847 0.892 -1.03 ; - comp large : flux_nf flux_f avecE0 sans E0 ; . phas'e 1.010 1.010 -0.467 ; . E0 seul phas'e 1.005 1.005 -0.450 ; . non phas'e 0.877 0.877 -0.450 ; => . c'est le dephasage de E0 qui produit un biais sur le flux, ; . retirer les branches ext. au pave ne change presque rien, ; . retirer la branhe int. de E0 (sur axe u) est catastrophique. ; Conclusions : ; . on ne filtre rien du tout, l'ecriture de FILTRE_PAVE est une con- ; nerie. On conserve plus bas l'ajustement de la comp continue ; sur l'image filtree mais on ne l'utilise pas. ; . le flux adopte est le flux par somme sur les bordures nulle, ; car le flux quadratique ne vaut que pour les sources compactes. flux_eq = total(im_int) ; sauvegarde du flux extrapolation quadratique. ; Calcul de la somme des bordures parallelement a l'axe v b_n = fltarr(npi) & b_s = fltarr(npi) ; bordures nord et sud. bande = fltarr(npi) ; pour essai for i = 0, npi-1 do begin b_n(i) = total (im_int (i, npi-1-(n_bord-1) : npi - 1)) b_s(i) = total (im_int (i, 0 : n_bord - 1)) bande(i) = total (im_int (i, npi/2 - n_bord/2 : npi/2 + n_bord/2)) endfor ; Controle des bordures icon = 0 ; a servi a trouver l'erreur d'inversion des indices. if ((i_stop eq 1) and (icon eq 1 )) then begin i_bande = 1 ; 1 pour tracer bande centrale. b_10 = bande / 10 if (i_bande eq 0) then begin ymax = max(b_n) > max(b_s) ymin = min(b_n) < min(b_s) endif if (i_bande eq 1) then begin ymax = max(b_n) > max(b_s) > max(b_10) ymin = min(b_n) < min(b_s) < min(b_10) endif dy = ymax - ymin y1 = ymin - 0.1 * dy & y2 = ymax + 0.1 * dy absc = indgen(npi) - npi/2 window, 0 plot, absc, b_n, xstyle=1, yrange=[y1, y2], ystyle=1, linestyle=2 oplot, absc, b_s - 0.005 * dy, linestyle=1 ; 0.005 pour separer. if (i_bande eq 1) then oplot, absc, b_10 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'bande central/10 ___, bordure N - -, bordure S ...' print, 'RH_MALC_IM_2D, avec n_bord =', n_bord print, ' c_c_p :', c_c_p tot_im = total(im_int) print, ' total (image) :', tot_im print, ' total de la bande centrale :', $ 10 * total(b_10) print, ' total des bordures N et S :', $ total(b_n), ' ', total(b_s) n_pts_bord = npi * n_bord ; nbre de pts dans une bordure moy = (total(b_n) + total(b_s)) / (2 * n_pts_bord) ; moy est la quantite dont il faut diminuer tous les points de ; l'image pour que les bordures aient une somme nulle. print, ' moyenne des bordures N et S :', moy corr_c_c = - npi*npi * moy ; correction de comp continue. print, ' correc de comp cont :', corr_c_c comp_cont = total(im_int - moy) ; total rectifi'e. ; comp_cont = total(im_int) + corr_c_c ; total rectifi'e. print, ' total corrige (par bandes nulles) :', comp_cont stop endif ; Calcul du flux par bordures nulles avant filtrage (pour comparaison). ; Le nombre de lignes et haut et en bas sur lesquelles on fait la moyen ; ne pour ajuster la base de l'image est definie en entree. ; Obsolete (avant avoir trouve l'erreur d'inversion des indices) : ; 10 est un compromis : ; - surestimamt de 8 % le flux d'une source compacte, ; - surestimant de 3% le flux d'une source large en (1-x^4)^4 de 60 ; pixels hors tout a sa base. moy = (total (im_int (*, 0 : n_bord-1)) + $ total (im_int (*, npi-1 - (n_bord-1) : npi-1)) ) / (2 * n_bord *npi) ; for i = 0, n_bord-1 do begin ; moy = (total (im_int(0:i, *)) + $ ; faux ; total (im_int(npi-1-i:npi-1, *))) / (2*n_bord*npi) ; faux. ; endfor ; Erreur corrigee le 7 mai 2002 : la boucle est inutile et les ; indices etaient invers'es. ; Ca donnait pourtant de bons resultats quand on utilisait E0 ; et ce n'etait pas par hasard : dire que la somme des bor- ; dures E et W devait etre nulle revenait a dire que l'image ; 1D (sans resolution en NS) devait tomber a zero au bord du ; champ, autrement dit que l'image 1D EW n'etait pas aliasee. ; c'est bien sur vrai avec E0 mais c'est faux sans E0, et ; dans ce dernier cas le resultat pouvait etre violemment ; faux, et d'une facon qui depend de la largeur de la source. im_int_nf = im_int - moy ; nf comme non filtre (sauvegarde avant le cal- ; cul ci-dessous sur l'image filtree im_int_2. flux_nf = total (im_int_nf) ; Rem : on obtiendra ci-dessous une meilleure estimation du flux et on ; ajoutera une constante a im_int_nf de facon que son nouveau total ; soit egal a cette meilleure valeur du flux. flux_nf_0 = flux_nf ; sauvegarde pour controle. ; Calcul d'un flux par integration sur un domaine ou l'image est notable ; Une evaluation imprecise de "moy" peut fortement affecter le flux ; obtenu par integration sur tout le champ interferomerique CI, ; puisque l'erreur sur le flux resulte de l'integration sur tout le ; CI, alors que le soleil est nettement moins etendu (et c'est encore ; plus vrai pour un centre d'orage tres compact). ; Pour reduire cet effet, on va integrer le flux de la source sur un ; domaine ou la source est notable (p ex. > 10% de son maximum), mais ; il faut evidemment ne pas exclure les secondaires negatifs proches ; du max d'une source compacte, ni les secondaires >0 oi <0 un peu ; lointains. ; On peut integrer : ; . tous les points de l'image sur un domaine limit'e autour du maxi- ; mum d'image, et dont on precrit l'extension, ; . sur tout le reste, la ou l'image en valeur absolue est > une frac- ; tion prescrite de son maximum. im_periph = im_int_nf ; sera mise a zero pres du max. max_im = max(im_int_nf, ip) jy = ip / npi ; indice de l'ordonnee du maximum. ix = ip mod npi ; --------- l'abscisse ---------- n_dom = nint(0.10 * npi) ; pour integrer sur une fraction du champ de ; part et d'autre du maximum d'image. ; cas du cygne : definit le domaine principal ; cas solaire : sera fortement complete par le ; domaine "peripherique" n1x = ix - n_dom & n2x = ix + n_dom n1y = jy - n_dom & n2y = jy + n_dom ; Precautions simples sur les indices (au 2 dec 2005 on ne gere pas la proxi- ; mite des bords du champ) if (n1x lt 0) then n1x = 0 & if (n2x gt npi - 1) then n2x = npi - 1 if (n1y lt 0) then n1y = 0 & if (n2y gt npi - 1) then n2y = npi - 1 ; Integration du flux proche im_proche = im_int_nf (n1x:n2x, n1y:n2y) flux_nf_proche = total(im_proche) ; Integration du flux lointain max_im = max(im_int_nf, ip) n_el_proche = n_elements(im_proche) ; pour controle seulement. im_periph(n1x:n2x, n1y:n2y) = 0 ; elimination du domaine du flux proche val_frac = 0.10 dom_periph = where(abs(im_periph) gt val_frac * max_im, n_el_loin) ; cas du cygne : contribution marginale au flux total ; cas solaire : conribution principale le plus souvent, compte tenu de ; la definition restreinte du domaine proche plus haut. if (n_el_loin eq 0) then flux_nf_periph = 0 if (n_el_loin gt 0) then flux_nf_periph = total(im_periph(dom_periph)) ; Calcul du flux total, restreint aux regions ou l'image est notable flux_nf_im_notable = flux_nf_proche + flux_nf_periph n_el_proche = n_elements(im_proche) ; Conclusion sur la determination du flux ; Essais le 2 dec 2005 sur : ; . cygne du 14 dec 2005 21:00 (sert a calibrer type II meme jour) ; . type II ----------- (dit "White") ; . soleil calme du 27 juin 2004 ; En prenant : ; . n_bord=10 sans RH_DPATCHFITS_NRH, ; . n_dom = nint(0.10 * npi) dans RH_MALC_IM_2D ; . val_frac = 0.10 ------------------ ; on trouve de determinations de flux (compact, flux_nf_im_notable etc ; qui semblent coherents a + ou _ 10%. En particulier on n'a plus de ; grosses erreur sur le flux "total" des sources compactes (qui pouvait ; etre <0 du fait du mauvais ajustement sur les bandes N et S). ; ; Cela permet d'ameliorer la base des images, en ajoutant a im_int_nf ; une constante telle que total(im_int_nf) = flux_nf_im_notable ; help, n_el_proche, n_el_loin ; help, flux_compact, flux_nf_proche, flux_nf_periph, flux_nf_im_notable ; help, flux_nf ; Amelioration de la base de l'image (d'apres la conclusion ci-dessus) correc_im = (flux_nf_im_notable - total(im_int_nf) )/ (npi^2) im_int_nf = im_int_nf + correc_im flux_nf = flux_nf_im_notable ; on ecrase le precedent flux_nf. La suite n'a pas besoin d'etre modi- ; fiee. On adopte un plus bas cette valeur pour flux_total (chercher ; "flux_total = flux_nf"). correc_im_rel = correc_im / max_im ; controle seulement. correc_flux_rel = (flux_nf - flux_nf_0) / flux_nf_0 ; controle seulement. if ( (version_solarsoft eq 0) and (i_stop eq 1) ) then begin ; 0 version test de Claude 1 SolarSoft ; permet de gerer impressions et stop print, ' RH_MALC_IM_2D :' ch_format = "(' amelioration base image : ', f5.2, ' %')" print, format=ch_format, 100 * correc_im_rel ch_format = "(' => ------------ estimation flux : ', f7.2, ' %')" print, format=ch_format, 100 * correc_flux_rel ; stop endif ; Controle image finale apres ajustement de la composante continue icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin i_sh_tv = 2 ; 1 shade_surf, 2 tvscl ; Trace du lobe if (i_sh_tv eq 1) then begin window, 2 shade_surf, congrid (im_int_nf, 512, 512, cubic=-0.5) endif if (i_sh_tv eq 2) then begin window, 2, xsize=512, ysize=512 tvscl, congrid (im_int_nf, 512, 512, cubic=-0.5) endif xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'Image apres ajustement de la composante continue' stop endif ; Fin controle image et lobe apres ajustement de la composante ; continue. ; DEBUT DE SEQUENCE INUTILE ; Calcul d'une image dont la TF est debarassee des branches debordant du ; pave central et (eventuellement) des harm EW de E0. Cette image ; (esperament moins pathologique que l'image initiale) est utilisee ; pour ajuster la composante continue. je0_ew = 1 ; 0 sans les harm de E0 interieurs au pave, 1 avec. ; C'est nettement plus mauvais avec je0_ew = 0. FILTRE_PAVE, im_int, je0_ew, $ ; entree im_int_f ; sortie ; Controle du filtrage de l'image icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin if (je0_ew eq 0) then ch_leg ='image filtree (sans harm de E0)' if (je0_ew eq 1) then ch_leg ='image filtree (avec harm de E0)' window, 0 shade_surf, im_int xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'image originale' window, 1 shade_surf, im_int_f xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_leg window, 2 shade_surf, im_int_f - im_int xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'image filtree - image initiale' stop endif ; SUITE DE SEQUENCE INUTILE ; Controle du lobe et se sa TF (pour rappel) icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin window, 0 tvscl, congrid(shift(abs(harl), npi/2, npi/2), 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(TF) du lobe theorique, centree' window, 1 tvscl, congrid(l_th, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Lobe theorique' window, 2 shade_surf, congrid(l_th, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Lobe theorique' stop endif ; SUITE DE SEQUENCE INUTILE ; Ajustement de la composante continue de l'image filtree for i = 0, n_bord-1 do begin ; i indice des lignes somm'ees. moy = (total(im_int_f(0:i, *)) + total(im_int_f(npi-1-i:npi-1, *))) / $ (2*n_bord*npi) endfor im_int_f = im_int_f - moy flux_f = total (im_int_f) ; On admets que le flux est le meme pour im_int ; et im_int_f. ; Essais ; im_int_3 = im_int - moy ; pour essais ; flux_f_3 = total (im_int_3) ; ; print, 'RH_MALC_IM_2D : flux_f =', flux_f, ' flux_f_3 =', $ ; flux_f ; stop ; => flux_f_3 est exactement egal a flux_f ; Controle icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_MALC_IM_2D : calcul du flux :' print, ' reduction de la TF au pave : oui non' print, ' extrapol quadr des bas harm Flux =', flux_eq print, ' total(bordures N et S) nul Flux =', flux_f, ' ', $ flux_nf print, ' (finalement adopte dans RH_MALC_IM_2D)' stop endif ; Remarque (16 mar 01). On appelle f_ex le flux par extrapolation ; de la TF et f_bn le flux deduit de la somme des bordures nulles. ; . Cas d'une source compacte : f_bn est moins precis que f_ex mais ; acceptable (qques %) pour un reseau phase, mais il devient vite ; tres faux (qques dixiemes) si E0 est dephasee. ; FIN DE SEQUENCE INUTILE ; Choix final du flux ; On adopte le flux calcul'e sans filtrage (voir plus haut conclusion ; des essais du 23 mars 01) im_int = im_int_nf flux = flux_nf flux_total = flux_nf harm(npi/2, npi/2) = flux_nf ; Il faut la comp continue reelle pour remplir le pave par interpola- ; tion autour de l'origine. On ecrase ici le flux compact. ; Rem du 11 dec 01 : le flux total, determine par la methode de la som- ; me de bordures N et S nulles, et beaucoup plus sensible aux ; erreurs de gain complexe des antennes que le flux compact, ob- ; tenu par extrapolation. Avec un Cygne simule et les gains ; "reguliers" le flux total est 0.66 et le flux compact 1.02 ; (au lieu de 1). ; Controle du flux et de l'image finale avec pave lacunaire (ou partiellement ; lacunaires avec les antennes AA apres le 26 dec 03). icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin ; Zoom sur la partie centrale des harmoniques (de -nh a + nh). nh = 4 h_i_c = abs(harm(npi/2-nh : npi/2+nh, npi/2-nh : npi/2+nh)) ; Trace du tableau des harmoniques et de son histogramme. amul = (freq / 164)^0.9 ; corriger spectre (cas du du Cygne). ; amul = 1 ; cas du Soleil. ; tvscl, congrid(amul * h_i_c, 512, 512, cubic=-0.5) ; window, 1 bin_size=0.01 ; plot, histogram(amul * h_i_c, binsize=bin_size), $ ; xrange=[0.1/bin_size, 1.5/bin_size] ; Calcul de flux ; amul = (freq / 164)^0.9 ; corriger spectre (cas du du Cygne). amul = 1 ; cas du Soleil. flux = total (im_int) flux_c = amul * flux ifreq = fix(freq) ch_format = "('Freq =', i4, ' MHz Flux =', f5.2," + $ "' Flux * freq / 164 =', f5.2)" print, format=ch_format, ifreq, flux, flux_c ; stop ; Trace de l'image i_sh_tv = 2 ; 1 shade_surf, 2 tvscl if (i_sh_tv eq 1) then begin window, 0 shade_surf, congrid (im_int, 512, 512, cubic=-0.5) endif if (i_sh_tv eq 2) then begin window, 0, xsize=512, ysize=512 tvscl, congrid (im_int, 512, 512, cubic=-0.5) endif xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'Image apres choix de l''image non filtree' stop endif ; Fin controle flux et image finale avec pave lacunaire. ; FIN DE TRAITEMENT D'IMAGE AVEC PAVE LACUNAIRE if (i_pave_plein eq 0) then goto, a3 ; vers fin de RH_MALC_IM_2D ; if (i_pave_plein gt 0) then print, 'RH_MALC_IM_2D : remplissage du pave' ; Detection du cas ou la polarisation est negative ; On risque d'avoir de problemes si la polar est negative : les phases ; au voisinage de l'origine sont 180 deg et non zero, et reparties de ; chaque cote => de possibles problemes d'interpolation si on inter- ; pole entre 2 determinations. ; On commence donc pa rendre l'image positive, on remplit le pave ; central, puis on change le signe de l'image et de sa TF avant de ; quitter RH_MALC_IM_2D. im_max = max(im_int, min = im_min) ; if (abs(im_max) gt abs(im_min)) then begin i_pos_neg = 1 ch_pos_neg = '(image laissee positive)' endif if (abs(im_max) lt abs(im_min)) then begin i_pos_neg = -1 ch_pos_neg = '(image rendue positive)' endif if (abs(im_max) eq abs(im_min)) then begin print, 'image avec meme valeur absolue pour mas et min' stop ; pour reflechir endif im_int = i_pos_neg * im_int ; sera recalculee avec le bon signe. harm = i_pos_neg * harm ; sera restaur'e apres remplissage pave if (i_pos_neg eq 0) then ch_chgt_signe = '' if (i_pos_neg eq -1) then ch_chgt_signe = ' changee de signe' ; DEBUT DE REMPLISSAGE DU PAVE LACUNAIRE (se fait apres le calcul de la cc) ; Recentrage de l'image sur son maximum pour obtenir une TF "quasi-reelle" et ; interpolable (au moins quand il n'y a pas plusieurs sources compactes ; sur le soleil) if (i_recentrage eq 0) then begin ; pour que ca soit defini en sortie meme dans ce cas. x_centr_auto = 0 y_centr_auto = 0 endif if (i_recentrage ge 1) then begin ; k_max et l_max sont specifies dans RH_DPNEW si i_recentrage=2. if (i_recentrage eq 1) then begin ; Recherche des coordonnees du maximum (qui seront directement en ca- ; naux interferometriques puisqu'on manipule l'image inter- ; ferometrique) if (ch_polar eq 'non polar') then begin max_im = max(abs(im_int), ip) k_max = (ip mod npi) ; ind de colonne (orig coin inf gauche) l_max = (ip / npi) ; ------ ligne k_max = k_max - npi/2 ; ind de colonne (origine au centre). l_max = l_max - npi/2 ; ------ ligne ----------------- k_max_auto = k_max ; pour passer au cas polar. l_max_auto = l_max ; ------------------------ endif if (ch_polar eq 'polar') then begin k_max_auto = k_max ; reconduction cas non polar ou choisi l_max_auto = l_max ; ----------------------------------- endif endif if (i_recentrage eq 2) then begin ; On calcule pas le decalage puisqu'on utilise le decalage ; choisi en 0.01 Rs dans RH_DPNEW sous les noms de x_centrage ; et y_centrag et convertis en canaux "interferometriques" ; dans RH_DPATCHFITS_NRH, ou on disposait de l'heure hmsc. endif print, 'RH_MALC_IM_2D : position de recentrage pr centre, ' + $ 'convertie en canaux interfe-' ch_format = "(18x, 'rometriques =', 2i5)" print, format=ch_format, k_max, l_max if (ch_polar eq 'polar') then begin print, ' (choisies ou reprises du cas non ' + $ 'polar automatique)' endif im_trans = shift(im_int, -k_max, -l_max) icon = 0 if (icon eq 1) then begin if (i_recentrage eq 1) then ch_auto = '(auto)' if (i_recentrage eq 2) then ch_auto = '(choisi)' ch_format = "(i5, i5, ' (canaux interfero).')" ch_recentrage = string(format=ch_format, k_max, l_max) im_512 = congrid(im_int, 512, 512, cubic=-0.5) im_t_512 = congrid(im_trans, 512, 512, cubic=-0.5) max_im = max(im_512) im_512 (256, *) = max_im & im_512 (*, 256) = max_im im_t_512 (256, *) = max_im & im_t_512 (*, 256) = max_im chsiz = 1.3 ; Image initiale window, 0, xsize=512, ysize=512 tvscl, im_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=chsiz, $ 'Image interferometrique' + ch_chgt_signe xyouts, 0.5, 0.93, /normal, alignment=0.5, charsize=chsiz, $ 'avant recentrage' ; Image recentree window, 1, xsize=512, ysize=512 tvscl, im_t_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=chsiz, $ 'Image interferometrique' + ch_chgt_signe xyouts, 0.5, 0.93, /normal, alignment=0.5, charsize=chsiz, $ 'recentrage ' + ch_auto + ' en : ' + ch_recentrage ; Abs(TF image) a_harm_512 = congrid(abs(harm), 512, 512) window, 2, xsize=512, ysize=512 tvscl, a_harm_512 xyouts, 0.5, 0.92, /normal, alignment=0.5, charsize=chsiz, $ 'abs (TF image)' stop wdel, 0 & wdel, 1 & wdel, 2 endif ; Calcul de la TF "quasi-reelle" (orig coin inf gauche) de l'image trans- ; latee. harm_1 = harm ; sauvegarde de la TF lacunaire non ; "quasi-reelle", centree en milieu ; de tableau. harm = shift(fft(shift(im_trans, -npi/2, -npi/2), 1), npi/2, npi/2) endif ; Sauvegarde du sous_pave central dans le cas des observations avec les anten- ; nes AA et de l'axe u, qui contiennent des harmoniques impairs vrais ; (observes) ; On laisse tomber les harm de rang 9 en NS, merdiques a gerer. h_axe_u = harm(*, npi/2) if (nbre_harm gt 576) then h_sous_pave = $ harm(npi/2 - 8 : npi/2 + 8, npi/2 - 8 : npi/2 + 8) icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin help, h_axe_u nz = 16 ; zooming de -nz a +nz. window, 0 absc = indgen(2*nz + 1) - nz plot, absc, abs(h_axe_u(npi/2 - nz : npi/2 + nz)), xstyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'abs(h_axe_u)' stop endif ; Controle de l'amplitude des harm observes et de ceux du lobe theorique icon = 0 ; changer icon aussi "apres remplissage". if ((icon eq 1) and (i_stop eq 1)) then begin print, 'MALC_IM_2D : amp des harmomiques avant remplissage du pave' nz = 6 ; zooming a_h_im = abs(harm(npi/2-nz:npi/2+nz, npi/2-nz:npi/2+nz)) a_h_l_th = abs(harl(npi/2-nz:npi/2+nz, npi/2-nz:npi/2+nz)) a_h_im_512 = congrid(a_h_im , 512, 512, cubic=-0.5) a_h_l_th_512 = congrid(a_h_l_th, 512, 512, cubic=-0.5) window, 0 shade_surf, a_h_l_th_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harmoniques du lobe theorique)' window, 1 shade_surf, a_h_im_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harmoniques observes)' window, 2 absc = indgen(2*nz) - nz plot, a_h_im(*, nz) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harmoniques observes sur axe u)' stop endif ; fin de controle ; Reconstitution des harm ew multiples impairs de 50m pour image et lobe ; Le remplissage du pave central se fait apres la reconstitution de la ; composante continue par la methode des bandes. Le but est en ; effet de desaliaser les images (eviter les distorsions par le "tun- ; nel, obtenir des flux par integration locale, etc.). ; Le remplissage va totalement desialiaser le lobe theorique, et le ; clean ne pourra donc plus desaliaser. Il faut donc que le remplis- ; sage le fasse, en particulier qu'il ne soit pas biaise pres de ; l'origine. ; Le remplissage du pave se faisant par interpolation, il faut que la ; cc soit correcte pour interpoler pres de l'origine. t0 = systime(1) ; temps en sec depuis le 1er jan 70 ; Interpolation bi-lineaire if (i_pave_plein eq 2) then begin ; interpolation lineaire. ; l'interpolation bi-lineaire correspondait a i_pave_plein=1 jusqu'au ; 9 oct 01. ; Calcul des tableaux des indices ix et iy des points du pave central ou on ; ne connait pas la visibilite ix = replicate(-1000, 33*45) iy = replicate(-1000, 33*45) ; On considere le pave amput'e de la 1ere et la derniere ligne, ; pour lesquelles l'interpolation ne peut se faire qu'en abscis- ; se. Le nombre de pts dans ce sous-pave est : 33*45. for j=0, 44 do begin for i=0, 30, 2 do begin ; 16 valeurs de i ix(i + 1 + j * 33) = -15 + i iy(i + 1 + j * 33) = -22 + j endfor endfor ; niveau de l'interpolation bi-lineaire. dom = where((ix ne -1000) and(iy ne 0)) ; exclut aussi l'axe x. ix = ix(dom) iy = iy(dom) icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin help, ix, iy, i, j stop print, ix stop print, iy stop endif ; verification OK ix = ix + npi/2 iy = iy + npi/2 ; Interpolation bilineaire pour l'image, remplissage avec 1 pour le lobe harm(ix, iy) = $ ( harm(ix-1, iy-1) + harm(ix-1, iy ) + harm(ix-1, iy+1) + $ harm(ix+1, iy-1) + harm(ix+1, iy ) + harm(ix+1, iy+1) ) / 6 harl(ix, iy) = $ ( harl(ix-1, iy-1) + harl(ix-1, iy ) + harl(ix-1, iy+1) + $ harl(ix+1, iy-1) + harl(ix+1, iy ) + harl(ix+1, iy+1) ) / 6 ; harl(ix, iy) = 1 for i=-15, 15, 2 do begin ; interpolation en x pour - et + 23. harm(npi/2+i, npi/2-23) = $ (harm(npi/2+i-1, npi/2-23) + harm(npi/2+i+1, npi/2-23)) / 2 harm(npi/2+i, npi/2+23) = $ (harm(npi/2+i-1, npi/2+23) + harm(npi/2+i+1, npi/2+23)) / 2 harl(npi/2+i, npi/2-23) = $ (harl(npi/2+i-1, npi/2-23) + harl(npi/2+i+1, npi/2-23)) / 2 harl(npi/2+i, npi/2+23) = $ (harl(npi/2+i-1, npi/2+23) + harl(npi/2+i+1, npi/2+23)) / 2 ; harl (npi/2+i, npi/2-23) = 1 ; harl (npi/2+i, npi/2+23) = 1 endfor endif ; fin de l'interpolation bi-lineaire. ; Interpolation cubic spline if (i_pave_plein eq 1) then begin ; interpolation cubic spline. ; l'interpolation cubic spline correspondait a i_pave_plein=2 jusqu'au ; 9 oct 01. ; Sauvegarde du sous-pave anti-aliasing dans le cas de la simulation if (i_sim_anti_alias eq 1) then begin i0 = 4 ; nbre d'antennes anti-aliasing. j0 = 8 ; rang de la NS ou est le sous-reseau EW ; i0 et j0 sont definis de facon identique dans SIMUL_RH. s_pave_0 = harm(np/2 - 2*i0 : np/2 + 2*i0, np/2 - j0 : np/2 + j0) endif ; Reduction a un pave reduit aux harmoniques pairs en EW, ie on saute les ; harmoniques impairs. ; L'interpolation avec CONGRID en doublera la dimension pour loger ; les harmoniques impairs) pave_lacun = complexarr(17, 47) for j=0, 46 do begin j_c = j - 23 ; indice centr'e dans harm. for i = 0, 32, 2 do begin i_c = i - 16 ; indice centr'e dans harm. pave_lacun ( 8 + i_c/2, 23 + j_c) = $ harm (npi/2 + i_c , npi/2 + j_c) endfor endfor ; Mesure du temps de calcul d'interpolation (on trouve 0.003 sec) t1 = systime(1) ch_1 = "('Temps de remplissage du pave " if (i_pave_plein eq 2) then ch_2 = "bilineaire'" if (i_pave_plein eq 1) then ch_2 = "spline'" ch_3 = "' sec')" ch_format = ch_1 + ch_2 + ", f6.3," + ch_3 ; print, format=ch_format, t1 - t0 ; Controle du pave reduit icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_MALC_IM_2D avant interpolation du pave reduit : ' print, ' abs(harm(npi/2, npi/2)) =', $ abs(harm(npi/2, npi/2)) print, ' abs(pave_lacun(16, 23)) =', $ abs(pave_lacun(16, 23)) window, 0 shade_surf, pave_lacun xyouts, 0.5, 0.97, /normal, alignment=0.5, 'abs(pave_lacun)' window, 1 tvscl, congrid(pave_lacun, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, 'abs(pave_lacun)' stop nz = 4 ; zooming pave_lacun de -nz a nz. absc_1 = (indgen(2*nz + 1) - nz) * 2 absc_2 = indgen(2*(2*nz) + 1) - 2*nz wset, 0 plot , absc_1, abs(pave_lacun(8-nz:8+nz, 23)), xstyle=1 oplot, absc_2, abs(harm(npi/2 -2*nz : npi/2 + 2*nz, npi/2)), $ linestyle=2 ; wset, 1 ; plot, absc_2, abs(harm(npi/2 -2*nz : npi/2 + 2*nz, npi/2)), $ ; xstyle=1 stop endif ; Calcul du pave complet'e par interpolation pave = congrid(pave_lacun, 34, 47, cubic=-0.99) ; cubic=-0.99 est le plus proche plus proche de l'interp theorique. ; Rem : la dimension en x doit etre exactement doublee (et non ; changee en 33 pour reconstituer les abscisses intermediaires. ; On laissera tomber les points les plus a droite, qui debordent ; de une rangee du champ du pave lacunaire (i=17). ; Rem : cubic=-0.5 oubli'e => interpolation d'ordre 0 (escalier) ; et image mauvaise. ; Rappel : si on a recentre l'image sur son maximum, alors a ce ; stade harm est quasi-reelle. pave = pave(0:32, *) if (i_sim_anti_alias eq 1) then begin s_pave_1 = pave(16 - 2*i0 : 16 + 2*i0, 23 - j0 : 23 + j0) ; pour evaluer la qualite de la reconsitution par interpolation endif ; Immersion harm_lacun = harm ; sauvegarde. harm(npi/2 - 16, npi/2 - 23) = pave ; immersion. axe_u_2 = harm(*, npi/2) ; pour controle. ; Controle de la qualite de l'interpolation (sur l'axe u) icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_MALC_IM_2D : controle interpolation sur axe u' nz = 8 ; zooming de -nz a +nz d_axe = axe_u_2 - h_axe_u absc = indgen(2*nz+1) - nz window, 0 & window, 1 wset, 0 plot , absc, abs(h_axe_u (npi/2-nz : npi/2+nz)), xstyle=1 oplot, absc, abs(axe_u_2 (npi/2-nz : npi/2+nz)), linestyle=2 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm(v=0) initial ___, interpole - -' wset, 1 plot, absc, abs(d_axe(npi/2-nz : npi/2+nz)) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm(v=0) interp[ole - (harm(v=0) initial' stop wset, 0 plot , absc, float(h_axe_u (npi/2-nz : npi/2+nz)), xstyle=1 oplot, absc, float(axe_u_2 (npi/2-nz : npi/2+nz)), linestyle=2 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'float(harm(v=0) initial ___, interpole - -' wset, 1 plot, absc, float(d_axe(npi/2-nz : npi/2+nz)) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'float(harm(v=0) interp[ole - (harm(v=0) initial' stop wset, 0 plot , absc, imaginary(h_axe_u (npi/2-nz : npi/2+nz)), xstyle=1 oplot, absc, imaginary(axe_u_2 (npi/2-nz : npi/2+nz)), linestyle=2 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'imaginary(harm(v=0) initial ___, interpole - -' wset, 1 plot, absc, imaginary(d_axe(npi/2-nz : npi/2+nz)) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'imaginary(harm(v=0) interp[ole - (harm(v=0) initial' stop endif ; Controle de l'allure de harm interpole hors de l'ave u ; Rappel : si on a recentre l'image sur son maximum, alors a ce ; stade harm est quasi-reelle. icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin i_ari = 1 ; traces : 0 abs, 1 float, 2 imaginary if (i_ari eq 0) then fonc = abs (harm) if (i_ari eq 1) then fonc = float (harm) if (i_ari eq 2) then fonc = imaginary(harm) window, 0 plot, absc, fonc(npi/2-nz : npi/2+nz, npi/2 + 1) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm(v=1))' window, 1 plot, absc, fonc(npi/2-nz : npi/2+nz, npi/2 + 2) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm(v=2))' window, 2 plot, absc, fonc(npi/2-nz : npi/2+nz, npi/2 + 3) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm(v=3))' window, 3 plot, absc, fonc(npi/2-nz : npi/2+nz, npi/2 + 4) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm(v=4))' stop endif ; fin de controle ; Restauration (eventuelle) du sous-pave anti_aliasing initial if (i_sim_anti_alias eq 1) then begin harm(npi/2 - 2*i0, npi/2 - j0) = s_pave_0 endif ; Comparaison sous-pave anti-aliasing initial et sous-pave interpole icon = 0 if((icon eq 1) and (i_stop eq 1) and (i_sim_anti_alias eq 1)) $ then begin print, 'RH_MALC_IM_2D : comparaison sous-pave avant et apres :' help, s_pave_0, s_pave_1 d_pav = s_pave_1 - s_pave_0 a_0 = abs(s_pave_0) a_1 = abs(s_pave_1) a_d = 100 * abs(d_pav) / max(abs(s_pave_0)) window, 0 shade_surf, a_0 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'abs(s_pave initial)' window, 1 shade_surf, a_1 xyouts, 0.5, 0.97, /normal, alignment=0.5, 'abs(s_pave interpole)' window, 2 shade_surf, a_d xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(s_pave ini - s_pave interp) , % de max( abs(s_pave ini))' stop endif ; Controle : visualisation de abs de harm et harm interpole icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_MALC_IM_2D : controle remplissage par cubic spline' nz = 12 ; zooming de-nz a + nz. ha_la = harm_lacun(npi/2-nz : npi/2+nz, npi/2-nz : npi/2+nz) ha = harm (npi/2-nz : npi/2+nz, npi/2-nz : npi/2+nz) ha_la_512 = congrid(ha_la, 512, 512, cubic=-0.5) ha_512 = congrid(ha , 512, 512, cubic=-0.5) window, 0 shade_surf, abs(ha_la_512) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm) avant completement' window, 1 shade_surf, abs(ha_512) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm) apres completement par interpolation' stop endif ; Complement du pave pour le lobe theorique for j = 0, 46 do begin j_c = j - 23 ; indice centr'e. for i = 0, 30, 2 do begin ; modif 15 harm impairs seuls. i_c = i - 15 ; indice centre harl (npi/2 + i_c, npi/2 + j_c) = 1 endfor endfor endif ; fin d'interpolation par cubic spline avec CONGRID. ; Restauration du sous-pave des harmoniques reellement observes avec les AA if (nbre_harm gt 576) then harm(npi/2 - 8, npi/2 - 8) = h_sous_pave ; Restauration des valeurs initiales de harm sur l'axe u si E0 est utilisee ; Rappel : si on a recentre l'image sur son maximum, alors a ce ; stade harm est quasi-reelle. if (ie0_ew eq 1) then harm(0, npi/2) = h_axe_u ; Essais : avec un soleil presque calme et un CME ; harm(*, npi/2) = 0 donne une gouttiere tres creuse. ; harm(*, npi/2) = h_axe_u --------- base ave 4.5 ondulations EW. ; harm(0, npi/2) = h_axe_u ----- strictement idem ; pas de restauration ----- base strictement plate ; Magouille sur l'harm 5, pour reduire la base ondulante quand le SC domine ; harm(npi/2 - 5, npi/2) = axe_u_2(npi/2 - 5) ; harm(npi/2 + 5, npi/2) = axe_u_2(npi/2 + 5) ; Controle de l'amplitude des harm observes et de ceux du lobe theorique icon = 0 ; changer icon aussi "avant remplissage". if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_MALC_IM_2D : amp des harmomiques apres remplissage du pave' nz = 6 ; zooming a_h_im = abs(harm(npi/2-nz:npi/2+nz, npi/2-nz:npi/2+nz)) a_h_l_th = abs(harl(npi/2-nz:npi/2+nz, npi/2-nz:npi/2+nz)) window, 2 shade_surf, a_h_l_th xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harmoniques du lobe theorique' window, 3 shade_surf, a_h_im xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harmoniques observes)' stop endif ; fin de controle ; Fin de remplissage du pave lacunaire sur la TF "quasi-reelle ; Restauration du signe initial de la TF de l'image ; (le signe a ete change pour un remplissage commode par interpolation ; du pave central) ; im_int = i_pos_neg * im_int ; inutile car recalculee plus bas harm = i_pos_neg * harm ; en commentaire pour essai seulement. ; Calcul de la TF (toujours centree en milieu de tableau) de l'image en sa ; position reelle. ; Cette TF n'est plus quasi-reelle. if(i_recentrage ge 1) then begin im_trans_plein = shift(float(fft(shift(harm, -npi/2, -npi/2), -1)), $ npi/2, npi/2) im_int = shift(im_trans_plein, k_max, l_max) ; Rappel : k_max et l_max sont les coordonnees de la position ; de recentrage par rapport au centre du champ. harm = shift(fft(shift(im_int, -npi/2, -npi/2), 1), npi/2, npi/2) ; cette valeur de harm ecrase la valeur quasi-reelle utilisee pour ; l'interpolation. endif ; Controle de l'image avec le pave plein icon = 0 if (icon eq 1) then begin ; le signe de l'image et de sa TF ayant ete restauree, on le change ; temporairement ici pour l'affichage de controle. window, 0, xsize=512, ysize=512 loadct, 0 tvscl, congrid (i_pos_neg * im_int, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'Image initiale, ' + ch_pos_neg window, 1, xsize=512, ysize=512 loadct, 0 tvscl, congrid (i_pos_neg * im_trans_plein, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'Image avec pave plein, ' + ch_pos_neg stop endif ; fin de controle ; Calcul d'image et de lobe theorique avec centre du Soleil en (npi/2, npi/2) if (i_recentrage eq 0) then begin im_int = float (shift(fft(shift(harm, -npi/2, -npi/2), -1), $ npi/2, npi/2)) ; inutile de refaire un calcul deja fait. endif l_th = float (shift(fft(shift(harl, -npi/2, -npi/2), -1), npi/2, npi/2)) ; Rem : total(im_int) = comp continue d'apres la propriete 2 des fft en ; IDL. Si la comp continue utilisee est mauvaise l'image est ; simplement decalee. flux = flux_nf ; reconduction du flux de l'image avec pave lacunaire. ; FIN DE CALCUL D'IMAGE AVEC PAVE PLEIN (la suite n'est que controles). ; Controle final de im_int et l_th (voir l'effet du remplissage du pave) icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_MALC_IM_2D : avant sortie' ; Trace lobe et im_int im_int_512 = congrid(im_int, 512, 512, cubic=-0.5) l_th_512 = congrid(l_th , 512, 512, cubic=-0.5) window, 0 shade_surf, l_th_512 if (i_pave_plein eq 0) then ch_2 = ', pave lacunaire' if (i_pave_plein gt 0) then ch_2 = ', pave plein' xyouts, 0.5, 0.976, /normal, alignment=0.5, $ 'l_th (sortie de malc_im_2d)' + ch_2 window, 1 shade_surf, im_int_512 xyouts, 0.5, 0.976, /normal, alignment=0.5, $ 'im_2d (sortie de malc_im_2d)' + ch_2 ; Coupes sur l'axe u window, 3 i_ari = 1 ; 0 abs, 1 float, 2 imaginary. if (i_ari eq 0) then begin fonc_0 = abs(h_axe_u) fonc_2 = abs(axe_u_2) ch_1 = 'abs' endif if (i_ari eq 1) then begin fonc_0 = float(h_axe_u) fonc_2 = float(axe_u_2) ch_1 = 'float' endif if (i_ari eq 2) then begin fonc_0 = imaginary(h_axe_u) fonc_2 = imaginary(axe_u_2) ch_1 = 'imagin' endif nz = 10 absc = indgen(2*nz) - nz plot , absc, fonc_0(npi/2 - nz : npi/2 + nz), xstyle=1 oplot, absc, fonc_2(npi/2 - nz : npi/2 + nz), linestyle=2 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_1 + $ '(harm sur axe u : observes __, interpoles - -)' stop ; Coupes en dehors de l'axe u wset, 0 fonc = abs(harm) plot , absc, abs(h_axe_u(npi/2-nz:npi/2+nz)) ; plot , absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 0) oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 1), linestyle=1 oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 2), linestyle=2 oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 3), linestyle=3 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm v : 0 __, 1 .., 2 - -, )' wset, 1 fonc = float(harm) plot , absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 0) oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 1), linestyle=1 oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 2), linestyle=2 oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 3), linestyle=3 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'float(harm v : 0 __, 1 .., 2 - -, )' window, 2 ; set, 2 fonc = imaginary(harm) plot , absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 0) oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 1), linestyle=1 oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 2), linestyle=2 oplot, absc, fonc(npi/2-nz:npi/2+nz, npi/2 + 3), linestyle=3 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'imagin(harm v : 0 __, 1 .., 2 - -, )' stop print, 'Fin de RH_MALC_IM_2D' endif ; fin du controle final de im_int et l_th. a3: J_M = 0 ; fin du saut de remplissage du pave. i_impression = 0 ; souvent utile en routine. if (i_impression eq 1) then begin print, 'RH_MALC_IM_2D sortie :' print, ' flux total =', flux_total , $ ' (par somme des bordures nulles)' print, ' flux compact =', flux_compact, $ ' (par extrapol quadratique de visibilite)' endif ; Controle d'image ajoute le 17 oct 2005 icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin i_sh_tv = 2 ; 1 shade_surf, 2 tvscl if (i_sh_tv eq 1) then begin window, 0 shade_surf, congrid (im_int, 512, 512, cubic=-0.5) endif if (i_sh_tv eq 2) then begin window, 0, xsize=512, ysize=512 tvscl, congrid (im_int, 512, 512, cubic=-0.5) endif xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'Image apres choix de l''image non filtree' stop endif ; Fin controle image finale. ; CONTROLES UTILES POUR LA CALIBRATION OU LA RECALIBRATION ; Trace des images 1D EW et NS ; utile controle de la re-calibration : il s'agit de tracer des images ; 1D de cygnes recalibres par eux-memes en utilisant les fichiers ; ceeff_corr* de facon exterieure a la procedure de calibration., ; ou de cygnes non recalibres. icon = 0 if (icon eq 1) then begin ch_freq = string(format="(i3, ' MHz : ')", nint(freq)) im_ew = total(im_int, 2) ; integration sur la 2eme dimension. im_ns = total(im_int, 1) ; ------------------ 1ere --------- nz = npi/4 ; zooming. im_ew_red = im_ew(npi/2-nz : npi/2+nz-1) im_ns_red = im_ns(npi/2-nz : npi/2+nz-1) im_ew_red_512 = congrid(im_ew_red, 512, 512, cubic=-0.5) im_ns_red_512 = congrid(im_ns_red, 512, 512, cubic=-0.5) fac_i = 512 / npi absc = (findgen(512) - 256) / fac_i * nz / (npi / 2) ; Trace image 1D EW window, 0, xsize=1000 amax = max(im_ew_red_512, min=amin) & da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da x1 = min(absc) - 2 & x2 = max(absc) + 2 plot, absc, im_ew_red_512, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_freq + $ 'image interferometrique 1D EW' ; Trace image 1D NS window, 1, xsize=1000 amax = max(im_ns_red_512, min=amin) & da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da x1 = min(absc) - 2 & x2 = max(absc) + 2 plot, absc, im_ns_red_512, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_freq + $ 'image interferometrique 1D NS' stop endif ; Trace des TF des images 1D EW et NS icon = 0 if (icon eq 1) then begin ch_freq = string(format="(i3, ' MHz : ')", nint(freq)) im_ew = total(im_int, 2) ; integration sur la 2eme dimension. im_ns = total(im_int, 1) ; ------------------ 1ere --------- tf_ew = shift(fft(shift(im_ew, -npi/2), 1), npi/2) tf_ns = shift(fft(shift(im_ns, -npi/2), 1), npi/2) amax = max(abs(tf_ew)) > max(abs(tf_ns)) tf_ew_1 = tf_ew + 0.0001 * amax tf_ns_1 = tf_ns + 0.0001 * amax a_ew = abs(tf_ew_1) a_ns = abs(tf_ns_1) p_ew = 180/!pi * atan(imaginary(tf_ew_1), float(tf_ew_1)) p_ns = 180/!pi * atan(imaginary(tf_ns_1), float(tf_ns_1)) ; zooming des amplitudes nz = npi/2 a_ew_red = a_ew(npi/2 : npi/2+nz-1) a_ns_red = a_ns(npi/2 : npi/2+nz-1) a_ew_red_512 = congrid(a_ew_red, 512, 512, cubic=-0.5) a_ns_red_512 = congrid(a_ns_red, 512, 512, cubic=-0.5) ; zooming des phases p_ew_red = p_ew(npi/2 : npi/2+nz-1) p_ns_red = p_ns(npi/2 : npi/2+nz-1) p_ew_red_512 = congrid(p_ew_red, 512, 512, cubic=-0.5) p_ns_red_512 = congrid(p_ns_red, 512, 512, cubic=-0.5) fac_i = 512 / (npi/2) absc = (findgen(512) ) / fac_i * nz / (npi / 2) ; Trace TF image 1D EW amplitude et phase window, 0, xsize=1000 amax = max(a_ew_red_512, min=amin) & da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da x1 = min(absc) - 2 & x2 = max(absc) + 2 plot, absc, a_ew_red_512, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_freq + $ 'abs(TF image EW)' window, 1, xsize=1000 amax = max(p_ew_red_512, min=amin) & da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da x1 = min(absc) - 2 & x2 = max(absc) + 2 plot, absc, p_ew_red_512, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_freq + $ 'phase(TF image EW)' stop ; Trace TF image 1D NS amplitude et phase window, 0, xsize=1000 amax = max(a_ns_red_512, min=amin) & da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da x1 = min(absc) - 2 & x2 = max(absc) + 2 plot, absc, a_ns_red_512, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_freq + $ 'abs(TF image 1D NS' window, 1, xsize=1000 amax = max(p_ns_red_512, min=amin) & da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da x1 = min(absc) - 2 & x2 = max(absc) + 2 plot, absc, p_ns_red_512, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_freq + $ 'phase(TF image 1D NS' stop endif ; Controle ultime d'image ajoute le 17 oct 2005 icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin i_sh_tv = 2 ; 1 shade_surf, 2 tvscl if (i_sh_tv eq 1) then begin window, 0 shade_surf, congrid (im_int, 512, 512, cubic=-0.5) endif if (i_sh_tv eq 2) then begin window, 0, xsize=512, ysize=512 tvscl, congrid (im_int, 512, 512, cubic=-0.5) endif xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'Image juste avant la sortie de RH_MALC_IM_2D' help, harl, l_th, harm, im_int, $ flux_total, flux_compact, uu, vv stop endif ; Fin controle image finale. ; Controle super-ultime ajoute le 21 fev 2006 : ; on visualise la tf reccalculee de l'image, apres la reconstitution ; de la composante continue. icon = 0 if ( (i_stop eq 1) and (icon eq 1) ) then begin taille = size(im_int) np2 = taille(2) / 2 tf_im = shift(fft(shift(im_int, -np2, -np2), 1), np2, np2) n_z = 32 n1z = np2 - n_z & n2z = np2 + n_z - 1 toto = tf_im(n1z : n2z, n1z : n2z) toto_512 = congrid(abs(toto), 512, 512) window, 0 loadct, 13 tvscl, toto_512 xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ 'abs(visibilite recalculee d''apres image) fin de RH_MALC_IM_2D' stop ; Le 21 fev 2006, on constate que la redondance est bien corrigee ; et qu'il n'y a pas ce creux selon l'axe v. ; On constate aussi que la composante continue a ete reconstituee. ; Ca a l'air OK. ; Le 23 fev 2006 on a trouve (apres une longue soiree d'efforts), ; une faute de cpier-clooer + fautre de frappe sur la correction ; des harmoniques de E0 dans CORR_GAIN_ANT, qui affectait les ; 8 1er harm NS de NS8. endif end ; fin de RH_MALC_IM_2D.