PRO IM_2D_CON, $ ; tous parametres en entree. sel_ant, $ i_stop, i_sys_lin, $ freq , ch_pol, ie0_ew, $ harm_N, harm_N_1, sref, correl, $ npi , n_bord, i_AA, $ icon_1, icon_2, icon_3, icon_4, $ n_win, frac_z, ch_siz, $ n_image, hmsc_deb, hmsc_fin ; Creation : 24 avril 03, a partir de IM_2D_CON ; But : calcul des images interferometrique 1D EW et NS pour les sources de ; controle. ; i_AA = 1 pour superposer les harm de NS8 avec les AA ; Modifications : ; 03 sep 26 : ajout de sel_ant a la liste d'entree. ; Il s'agit, dans le cas polar ou les harmoniques de NS8 sont ; nuls (sauf avec NS0_ew), de compenser la correction de redon ; dance, sui divise par 3. ; Je prefere mettre la rustine ici, ou sont seules concernees les ; observations de controle, plutot que RH_HARM_N_VIS, qui traite ; aussi le soleil. ; 04 jan 15 extansion aux observations avec lea antennes AA (apres jan 04). ; 05 nov 24 : n_image, hmsc_deb, hmsc_fin ajoutes pour affichage sur le ; trace ; ETAPES IM_2D_CON ; - 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). ; - centrage (dans le cas du soleil) pour avoir des TF quasi_relles. ; - les images 1D EW et NS sont obtenus par integration de l'image 2D sur ; les directions NS et EW respectivement. ; Notations: ; nbre_harm nombre d'harmoniques (576 sans les AA, 648 avec les AA). ; nbre_ant nombre d'antennes (44 ou 48). ; 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_N_1 sauvegarde de l'harmonique (NS0_ns, E1) avant mise a zero ; eventuelle. ; harm, harl complexarr (128, 128) tableaux des harm dans le plan (u, v). ; En sortie ils sont sentres 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. ; sref 'soleil' ou 'autre' ; uu, vv coord de frequences spatiales dans le plan (u, v). ; print, 'IM_2D_CON entree : k_max, l_max =', k_max, l_max nbre_harm = n_elements(harm_N) if (nbre_harm eq 576) then nbre_ant = 44 if (nbre_harm gt 576) then nbre_ant = 48 ; 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='IM_2D_CON : 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, 'IM_2D_CON : entree :' ; Controle des harm NS avec NS0 icon = 0 if (icon eq 1) then begin harm_ns = complexarr(25) for j = 1, 24 do harm_ns(j) = harm_N(9 + 18*j) nz = 24 ; zoom de 0 a nz. NUMERO D'ANTENNE ET NON D'HARMONIQUE window, 0 plot, abs(harm_ns(0:nz)) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm avec NS0_ew)' stop endif ; REMPLISSAGE DU PLAN UV AVEC LES HARMONIQUES OBSERVES i_stop = 0 ; i_stop introduit dans rh_harm_n_vis le 22 fev 06. RH_HARM_N_VIS, $ ; fichier RH_HARM_N_VIS.PRO. i_stop, $ ; entree npi, correl, harm_N, $ ; entree harm, harl, uu, vv ; sortie ; Produit les repartitions "harm" et "harl" des harmoniques spatiaux ; 2D, corrigees des redondances et CENTRES EN MILIEU DE TABLEAU, ; respectivement 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. Pour le Cygne polar ; les harm NS purs de 1 a 8 sont alors sous-evalues si NS8 est uti- ; lisee, puisque ses harmoniques polar sont nuls sauf avec NS0_ew, ; voir gestion de c cas plus bas. ; Gestion du cas polar : la correction de redondance dans RH_HARM_N_VIS ; a divis'e par 3 les harmoniques NS purs de 1 a 8 si NS8 est utilisee. ; Or les harm NS de NS8 sont nuls en polar sur le Cygne, sauf avec ; NS0_ew), et tous ces harmoniques ont ete divises par 3. Il faut : ; . multiplier par 3 les harmoniques de 1 a 7 (dans la redondance ; triple seuls les harm de NS0_ew avec les NS sont non nuls ; . multiplier par 2 l'harm 8 puisque (NS8, NS16) est nul ; . multiplier par 2 l'harm 9 puisque (NS8, NS17) est nul ; Bien sur cela n'arrange que les traces 1D NS. n_NS8 = 27 ; numero de NS8 if ((ch_pol eq 'polar') and (sel_ant(n_NS8) eq 1)) then begin print, 'IM_2D_CON : compensation de la redondance pour les harm NS' harm(npi/2, npi/2 - 7) = harm(npi/2, npi/2 - 7 : npi/2 + 7) * 3 harm(npi/2, npi/2 - 8) = harm(npi/2, npi/2 - 8) * 2 harm(npi/2, npi/2 + 8) = harm(npi/2, npi/2 + 8) * 2 harm(npi/2, npi/2 - 9) = harm(npi/2, npi/2 - 9) * 2 harm(npi/2, npi/2 + 9) = harm(npi/2, npi/2 + 9) * 2 endif ; Controle amplitude de harm et harl icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'IM_2D_CON : 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 avant les ; AA et je ne me fais pas chier a modifier pour le cas ou on les a). 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 pricipe meilleure evaluation du flux ; compact peut quelqufois planter. ; Avait ete mis a zero dans RH_MALC_IM_2D sur ; la requete de D. Jacquin pour eviter ; plantages dans traitement s de masse. if (i_extrapol_quad eq 0) then begin print, 'IM_2D_CON : pas d''extrapolation quadratique du flux' endif 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 IM_2D_CON 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 IM_2D_CON' 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. ; 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. ; 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 : ; - bew et bns interviennent avec un signe + dans phi (voir plus haut, ; modif du 27 mai 99). ; - l'amplitude de la composante continue et des bas harmoniques est ; egale au flux (sfu). ; Conclusions pour le calcul de l'image : ; - la TF inverse (drapeau 1) realise l'operation faite par le recep- ; teur du RH, avec une composante continue egale au flux en sfu ; (propriete 1) et des bas harmoniques qui en sont voisins. ; - la TF directe (drapeau -1) calcule l'image de la source, normalisee ; de sorte que la somme des points soit egale a son flux (sfu) a ; partir des harmoniques fournis par le RH (propriete 2). ; En pratique on ne calcule le flux que quand l'image est exprimee ; en sfu (et non en Kelvins). ; Calcul d'image et de lobe theorique avec origine 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. ; Fin de calcul d'image avec pave lacunaire ; 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, 'IM_2D_CON, 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) 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) ; ne sert qu'a un controle. ; Choix final du flux ; Pour les sources de controle le flux par bordures n'est pas fiable ; car il est sensible aux erreurs de gain complexe des antennes. ; On adopte le flux compact (parabolique). im_int = im_int_nf flux = flux_nf flux_total = flux_nf harm(npi/2, npi/2) = flux_compact ; Controle du flux et de l'image finale avec pave lacunaire. 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 & wset, 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 window, 2 & wset, 2 tvscl , congrid (im_int, 512, 512, cubic=-0.5) shade_surf, congrid (im_int, 512, 512, cubic=-0.5) stop window, 0 & wset, 0 ; pour la suite. stop endif ; Fin controle flux et image finale avec pave lacunaire. i_impression = 0 ; souvent utile en routine. if (i_impression eq 1) then begin print, 'IM_2D_CON sortie :' print, ' flux total =', flux_total , $ ' (par somme des bordures nulles)' print, ' flux compact =', flux_compact, $ ' (par extrapol quadratique de visibilite)' endif ; Recentrage de la source dans le cas du soleil pour avoir des TF quasi-reelles ; On recentre sur le max de l'image 2D if (sref eq 'SOLEIL') then begin np2 = npi ; si np2 > npi la phase des harmoniques impairs au dela de 32 ;n'est plus nulle. im_np2 = congrid(im_int, np2, np2, cubic=-0.5) ; interpolation image pour avoir precision sur position maximum. max_im = max(abs(im_np2), ip) k_max = (ip mod np2) ; ind de colonne (orig coin inf gauche) l_max = (ip / np2) ; ------ ligne kmax = k_max - np2/2 ; abscisse du max par rapport au centre lmax = l_max - np2/2 ; ordonnee ---------------------------- im_np2 = shift(im_np2, -kmax, -lmax) im_int = congrid(im_np2, npi, npi, cubic=-0.5) endif ; fin de recentrage dans le cas du soleil. ; Definition de chaines de caracteres pour affichage sur les traces ch_freq = string(format="(i3, ' MHz : ')", nint(freq)) ch_pol_freq = ch_pol + ' ' + ch_freq ; definition de ch_integ ch_n_image = string(n_image) if (n_image eq 1) then ch_image = ' image, de ' if (n_image gt 1) then ch_image = ' images, de ' ch_format = "(i2.2, ':', i2.2, ':', i2.2)" ch_integ = 'integration sur ' + strcompress(string(n_image), $ /remove_all) + ch_image + $ string(format=ch_format, hmsc_deb(0:2)) + ' a ' + $ string(format=ch_format, hmsc_fin(0:2)) ; Detection d'une image avec max absolu negatif (possible en polar) max_im = max(im_int, min=min_im) if (abs(max_im) gt abs(min_im)) then i_pos_neg = 1 if (abs(max_im) lt abs(min_im)) then i_pos_neg = -1 if (abs(max_im) eq abs(min_im)) then begin print, 'image avec max et min egaux en valeur absolue' stop ; reflechir. endif if (i_pos_neg eq 1) then begin ch_chgt_signe = 'non changee de signe,' endif if (i_pos_neg eq -1) then begin harm = - harm im_int = - im_int ch_chgt_signe = 'changee de signe,' endif ; Trace des images 1D EW et NS des sources de controle ; utile controle de la re-calibration : il s'agit de tracer des images ; 1D de Cygnes recalibres par eux-memes en utilisant les fichiers ; coeff_corr* de facon exterieure a la procedure de calibration, ; ou de Cygnes non recalibres. if (icon_1 eq 1) then begin im_ew = total(im_int, 2) ; integration sur la 2eme dimension. im_ns = total(im_int, 1) ; ------------------ 1ere --------- nz = nint(0.5 * npi) ; 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) ; Dilatation par un facteur 10 im_ew_cen = im_ew(npi/2 - nz : npi/2 + nz - 1) ; interpolation a 512 im_ew_red_512 = congrid(im_ew_red, 512, cubic=-0.5) im_ns_red_512 = congrid(im_ns_red, 512, cubic=-0.5) ; Dilatation par un facteur 10 im_ew_cen = congrid(im_ew_red_512, 5120, cubic=-0.5) im_ns_cen = congrid(im_ns_red_512, 5120, cubic=-0.5) im_ew_cen_512 = im_ew_cen(2560 - 256 : 2560 + 256 - 1) im_ns_cen_512 = im_ns_cen(2560 - 256 : 2560 + 256 - 1) fac_i = 512 / npi absc = (findgen(512) - 256) / fac_i * nz / (npi / 2) ; Trace image 1D EW loadct, 0 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 oplot, absc, im_ew_cen_512, linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=ch_siz, $ ch_freq + 'image interferometrique 1D EW. Dilatation par 10 : ...' xyouts, 0.5, 0.93, /normal, alignment=0.5, charsize=ch_siz, $ ch_chgt_signe + ' ' + ch_integ ; Trace image 1D NS loadct, 0 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 oplot, absc, im_ns_cen_512, linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=ch_siz, $ ch_freq + 'image interferometrique 1D NS. Dilatation par 10 : ...' xyouts, 0.5, 0.93, /normal, alignment=0.5, charsize=ch_siz, $ ch_chgt_signe + ' ' + ch_integ stop absc = 0 ; precaution pour la suite endif ; Preliminaire a l'ajout de l'harm (E1, NS0_ns) au trace de TF (image EW) a_E1_NS0_ns = abs(harm_N_1) p_E1_NS0_ns = -180/!pi * atan( imaginary(harm_N_1), float(harm_N_1) ) ; signe - pour considerer E1 comme reference. ; Preliminaire a l'ajout des harmoniques de NS8 avec les AA (bases E -> W) if (nbre_harm gt 576) then begin h_AA1 = harm_N(576 + 8) & h_AA2 = harm_N(594 + 8) h_AA3 = harm_N(612 + 8) & h_AA4 = harm_N(630 + 8) a_h_AA1 = abs(h_AA1) & a_h_AA2 = abs(h_AA2) a_h_AA3 = abs(h_AA3) & a_h_AA4 = abs(h_AA4) p_h_AA1 = 180/!pi * atan(imaginary(h_AA1), float(h_AA1) + 0.00001) p_h_AA2 = 180/!pi * atan(imaginary(h_AA2), float(h_AA2) + 0.00001) p_h_AA3 = 180/!pi * atan(imaginary(h_AA3), float(h_AA3) + 0.00001) p_h_AA4 = 180/!pi * atan(imaginary(h_AA4), float(h_AA4) + 0.00001) endif ; Trace des TF des images 1D EW et NS 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)) ; recadrage rustineux des phases dom = where(p_ew gt 100, ic) if (ic gt 0) then p_ew(dom) = p_ew(dom) - 360 ; zooming des amplitudes ; utile : zoom differents en EW et NS. nz_ew = 62 ; usuel npi/2 nz_ns = 45 ; usuel npi/2 a_ew_red = a_ew(npi/2 : npi/2 + nz_ew - 1) a_ns_red = a_ns(npi/2 : npi/2 + nz_ns - 1) ; zooming des phases p_ew_red = p_ew(npi/2 : npi/2 + nz_ew - 1) p_ns_red = p_ns(npi/2 : npi/2 + nz_ns - 1) ; Exclusion de la composante continue, qui est une invention de IM_2D_CON ; on exclut le 1er point non interpole, cad les 8 1ers points ; interpoles a 512. absc_ew = indgen(nz_ew - 1) + 1 ; sans interpoler a 512. absc_ns = indgen(nz_ns - 1) + 1 ; --------------------- b_ew_red = a_ew_red (1 : nz_ew - 1) b_ns_red = a_ns_red (1 : nz_ns - 1) q_ew_red = p_ew_red (1 : nz_ew - 1) q_ns_red = p_ns_red (1 : nz_ns - 1) a_r = 1 ; coeff de reduction de l'echelle de trace des amplitu- ; des, pour que l'ehelle ne soit pas imposee par les ; bas harmoniques dans le cas ou le soleil calme domi- ; ne. ; usuel : 1 pour les sources, 5 ou 10 pour le soleil. ; Trace TF image 1D EW amplitude et phase if (icon_2 eq 1) then begin window, 0, xsize=1000 amax = max(b_ew_red) > a_E1_NS0_ns amin = min(b_ew_red) < a_E1_NS0_ns if ((i_AA eq 1) and (nbre_harm gt 576)) then begin amax = amax > a_h_AA1 > a_h_AA2 > a_h_AA3 > a_h_AA4 amin = amin < a_h_AA1 < a_h_AA2 < a_h_AA3 < a_h_AA4 endif da = amax - amin y1 = amin - 0.05 * da & y2 = (amax + 0.05 * da) / a_r x1 = min(absc_ew) - 2 & x2 = max(absc_ew) + 2 plot, absc_ew, b_ew_red, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 oplot, [16], [a_E1_NS0_ns], psym=2 if ((i_AA eq 1) and (nbre_harm gt 576)) then $ oplot, [1, 3, 5, 7], [a_h_AA1, a_h_AA2, a_h_AA3, a_h_AA4], psym=4 ; psym=2 *, 3 ., 4 losange xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=ch_siz, $ ch_freq + 'abs de : TF(image EW ' + ch_chgt_signe + $ ') __ , abs(E1, NS0_ns) *, abs(NS8, AA) losange ' window, 1, xsize=1000 amax = max(q_ew_red) > p_E1_NS0_ns amin = min(q_ew_red) < p_E1_NS0_ns if ((i_AA eq 1) and (nbre_harm gt 576)) then begin amax = amax > p_h_AA1 > p_h_AA2 > p_h_AA3 > p_h_AA4 amin = amin < p_h_AA1 < p_h_AA2 < p_h_AA3 < p_h_AA4 endif da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da x1 = min(absc_ew) - 2 & x2 = max(absc_ew) + 2 plot, absc_ew, q_ew_red, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 oplot, [16], [p_E1_NS0_ns], psym=2 if ((i_AA eq 1) and (nbre_harm gt 576)) then $ oplot, [1, 3, 5, 7], [p_h_AA1, p_h_AA2, p_h_AA3, p_h_AA4], psym=4 ; psym=2 *, 3 ., 4 losange xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=ch_siz, $ ch_freq + 'phases de : TF(image EW ' + ch_chgt_signe + $ ') ___ , (E1(ref), NS0_ns) *, (NS8, AA) losange' stop endif ; Trace TF image 1D NS amplitude et phase loadct, 0 ; 0 BW, 3 red temp., 13 rainbow if (icon_3 eq 1) then begin window, 0, xsize=1000 amax = max(b_ns_red, min=amin) & da = amax - amin y1 = amin - 0.05 * da & y2 = (amax + 0.05 * da) / a_r x1 = min(absc_ns) - 2 & x2 = max(absc_ns) + 2 plot, absc_ns, b_ns_red, xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=ch_siz, $ ch_freq + 'abs(TF image 1D NS ' + ch_chgt_signe + ')' window, 1, xsize=1000 ; Recadrage rustineux pour (parfois mieux voir le saut de phase de NS45) dom = where(q_ns_red gt 120, ic) if (ic gt 0) then q_ns_red(dom) = q_ns_red(dom) - 360 ; Fin de recadrage rustineux amax = max(q_ns_red, min=amin) & da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da x1 = min(absc_ns) - 2 & x2 = max(absc_ns) + 2 plot, absc_ns, q_ns_red , xrange=[x1, x2], xstyle=1, $ yrange=[y1, y2], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=ch_siz, $ ch_freq + 'phase(TF image 1D NS ' + ch_chgt_signe + ')' stop wdel, 0 & wdel, 1 endif ; Trace image 2D if (icon_4 eq 1) then begin i_stop = 0 i_pave_plein = 0 i_recentrage = 0 ; recentrage d'image avant interpolation de TF, ; 0 sans recentrage ; 1 avec ---------- sur son maximum d'image ; 2 ----------------sur coord interfero ; ci-dessous. x_centrage = 0 ; absc interferometrique en canaux p r centre ; image. y_centrage = 0 ; ord -------------------------------------- RH_MALC_IM_2D, $ ; entree: i_stop, i_sys_lin, $ ch_pol, $ i_pave_plein, i_recentrage, x_centrage, y_centrage, $ freq, ie0_ew, harm_N, correl, npi, n_bord, $ ; sortie: image et lobe theorique, centre Soleil en (npi/2, npi/2). x_centrage_np, y_centrage_np, $ tf_l_th, l_th, tf_im_2d, im_2d, $ flux_total, flux_compact, uu, vv ; uu, vv coord de freq spatiales dans le plan (u, v), ; ne servent qu'a des delouiseries. ; tf_l_th, tf_im_2d sont centrees en milieu de tableau. ; ajoutes pour compatibilite le 13 mars 03. im_2d = i_pos_neg * im_2d ; Verification des flux totaux et compacts icon_5 = 0 if (icon_5 eq 1) then begin print, 'RH_CALC_IMAGE_HELIO : sortie de RH_MALC_IM_2D' print, ' flux total =', flux_total, $ ; ' (par somme des bordures nulle)' print, ' flux compact =', flux_compact, $ ; ' (par extrapol quadratique de visibilite)' stop endif ; Trace de im_2d et l_th (voir l'effet du remplissage du pave) ch_format = "('IM_2D_CON : sortie de RH_MALC_IM_2D, max d " + $ "image =', f7.3)" print, format=ch_format, max(im_2d) ; Zooming frac_c = frac_z(0) ; facteur de zooming sur le champ nz = nint(frac_c * npi/2) & n1 = npi/2 - nz & n2 = npi/2 + nz - 1 im_2d_512 = congrid(im_2d(n1:n2, n1:n2), 512, 512, cubic=-0.5) l_th_512 = congrid(l_th (n1:n2, n1:n2), 512, 512, cubic=-0.5) if (i_pave_plein eq 0) then ch_2 = 'pave lacunaire, ' if (i_pave_plein gt 0) then ch_2 = 'pave plein, ' ; Choix de l'echelle de couleur pour tvscl j_couleur = 13 ; 0 BW, 3 red temp., 13 rainbow ; Lobe theorique window, n_win, xsize=512, ysize=512 loadct, j_couleur tvscl , l_th_512 loadct, 0 xyouts, 0.5, 0.976, /normal, alignment=0.5, charsize=ch_siz, $ ch_pol_freq + 'l_th (sortie de malc_im_2d)' + ch_2 ; Image (en pratique image du Cygne, plus rarement du soleil) window, n_win + 1, xsize=512, ysize=512 loadct, j_couleur tvscl, im_2d_512 loadct, 0 xyouts, 0.5, 0.976, /normal, alignment=0.5, charsize=ch_siz, $ ch_pol_freq + 'image interferometrique ' + ch_chgt_signe xyouts, 0.5, 0.93, /normal, alignment=0.5, charsize=ch_siz, $ ch_2 + ch_integ i_saut = 0 ; 1 pour sauter shade_surf. if (i_saut eq 0) then begin azim = 30 ; rotation shade_surf (degres sens horaire). window, n_win + 2 loadct, 0 ; 0 BW, 3 red temp., 13 rainbow shade_surf, im_2d_512, az=azim loadct, 0 xyouts, 0.5, 0.976, /normal, alignment=0.5, charsize=ch_siz, $ ch_pol_freq + 'image interferometrique ' + ch_chgt_signe xyouts, 0.5, 0.93, /normal, alignment=0.5, charsize=ch_siz, $ ch_2 + ch_integ ; xyouts, 0.5, 0.976, /normal, alignment=0.5, charsize=ch_siz, $ ; ch_pol_freq + 'image interferometrique' + ch_2 ; xyouts, 0.5, 0.93, /normal, alignment=0.5, charsize=ch_siz, $ ; ch_integ endif stop ; wdel, n_win + 1 & wdel, n_win + 4 ; Trace de abs(TF) de l'image (voir si on utilise NS45) tf_im = shift(float(fft(shift(im_2D, -npi/2, -npi/2), 1)), npi/2,npi/2) ; Zooming frac_uv = frac_z(1) ; facteur de zooming dans le plan uv. tf_im_0 = tf_im ; sauvegarde if (frac_uv lt 0.95) then begin nc = npi/2 ; comme "centre" nz = nint(frac_uv * npi/2) ; comme "zoom" n1 = nc - nz & n2 = nc + nz - 1 tf_im = tf_im(n1 : n2, n1 : n2) endif a_tf = abs(tf_im) toto = tf_im + 0.0001 * max(a_tf) ; precaution pout atan. p_tf = 180/!pi * atan(imaginary(toto), float(toto)) ; a_tf_512 = congrid(a_tf, 512, 512, cubic=-0.5) a_tf_512 = congrid(a_tf, 512, 512) ; Visualisation de la croix fac_i = 512 / (128 * frac_uv) ; facteur d'interpolation de TF zoomee. n_trans = nint(fac_i / 2) ; n_trans = 0 ; cadavre laisse a tout hasard. a_tf_512_croix = a_tf_512 a_tf_512_croix (*, 256 + n_trans) = max(a_tf) a_tf_512_croix (256 + n_trans, *) = max(a_tf) ; pour vraiment centrer la croix. ; Visualisation des limites du pave n_u_1 = 256 + n_trans - nint(fac_i * 16) n_u_2 = 256 + n_trans + nint(fac_i * 16) n_v_1 = 256 + n_trans - nint(fac_i * 23) n_v_2 = 256 + n_trans + nint(fac_i * 23) if ((n_v_1 ge 0) and (n_v_2 le 512)) then begin a_tf_512_croix (n_u_1 : n_u_2, n_v_1) = max(a_tf) a_tf_512_croix (n_u_1 : n_u_2, n_v_2) = max(a_tf) a_tf_512_croix (n_u_1, n_v_1 : n_v_2) = max(a_tf) a_tf_512_croix (n_u_2, n_v_1 : n_v_2) = max(a_tf) endif ; Mise a zero de l'origine pour bien la voir sur le trace (idiot). ; a_tf_512(256, 256) = 0 & a_tf_512(257, 256) = 0 ; a_tf_512(257, 257) = 0 & a_tf_512(256, 257) = 0 ; a_tf_512(255, 257) = 0 & a_tf_512(255, 256) = 0 ; a_tf_512(255, 255) = 0 & a_tf_512(256, 255) = 0 ; a_tf_512(257, 255) = 0 loadct, 0 ; 0 BW, 13 rainbow A = 2000 toto = alog(1 + A * a_tf_512_croix / max (a_tf_512_croix)) ch_A = string(format="('A =', i2)", A) window, n_win + 4 loadct, 0 ; 0 BW, 3 red temp., 13 rainbow shade_surf, a_tf_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ ch_pol_freq + ' : abs(TF) de l''image' window, n_win + 5, xsize=512, ysize=512 loadct, 0 ; 0 BW, 3 red temp., 13 rainbow tvscl, toto xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ ch_pol_freq + ' : abs(TF) de l''image' xyouts, 0.5, 0.94, /normal, alignment=0.5, charsize=1.3, $ ' saturation en log(1 + A x/xmax) avec ' + ch_A stop endif ; fin du if sur icon_4 end ; fin de IM_2D_CON. ;______________________________________________________________________________ PRO RH_IMAGE_CON, $ ; tous les parametres sont en entree. ; Caracterisation de l'image : sref, correl, freq, ch_pol, npi, n_bord, $ ; Suppression (eventuelle) des harmoniques de certaines antennes ie0_ew, ie0_ns, ie1_ew, ie2_ew, ie2_ns, $ ins45_ew, ins45_ns, i_redond, gap_ns0, sel_ant, $ ; Donnees pour une image : n_image, hmsc_deb, hmsc_fin, $ harm_N, $ ; Divers i_sys_lin, i_AA, $ ; Affichages dans IM_2D_CON icon_1, icon_2, icon_3, icon_4, n_win, frac_z, ch_siz ; Creation 24 avril 03, a partir de RH_CALC_IMAGE_HELIO et RH_MALC_IM_2D. ; Fonction : tracer des images 1D des sources de controle ; - suppression eventuelle de certains harmoniques, ; - controle de la distribution des harmoniques dans "harm_N", ; - calcul de l'image (non polar ou polar) "sale" sur la grille des coord ; interferometriques entieres et normalisation, ; Modifications : ; 03 oct 10: - extension a l'utilisation des antennes AA. ; 05 nov 24: - n_image, hmsc_deb, hmsc_fin sont ajoutes pour affichage ; sur le trace. ; i_AA = 1 pour superposer les harm de NS8 avec les AA. ; Notations: ; harm_N vecteur des correlations. Rempli en non-polar ou polar juste ; avant l'appel a RH_IMAGE_CON. ; nbre_harm nombre de correlations (576, 648 ou 720) dans harm_N. ; correl tableau (2, nbre_harm) indiquant les numeros des 2 antennes corres- ; pondant a chacun des nbre_harm harmoniques. Permet, dans ; RH_MALC_IM_2D le passage de harm_N au tableau "harm" des harmo- ; niques a 2 dimensions. ; ; npi dimension de l'image interferometrique, figee a 128 apres le 6 jul ; 01. ; sref 'soleil' ou 'autre' ; Discernement des observations avec lea antennes anti-alias nbre_harm = n_elements(harm_N) ; Controle repris de RH_CALIB_CON (pour verif calib polar) icon = 0 if (icon eq 1) then begin print, 'RH_IMAGE_CON entree' hNS_NS0_ew = complexarr(25) for i=1, 24 do hNS_NS0_ew(i) = harm_N(9 + 18 *i) ch_format ="(50f8.3)" print, ' ' print, 'NS0_ew avec les NS' print, format=ch_format, abs(hNS_NS0_ew) print, ' ' print, 'NS0_ns avec NS7 a NS23' print, format=ch_format, abs(harm_N(573:575)), abs(harm_N(558:572)) print, ' ' print, 'NS8 avec NS0 a NS17' print, format=ch_format, abs(harm_N(540:557)) print, ' ' print, 'NS45 avec les NS' print, format=ch_format, abs(harm_N(450:485)) print, ' ' stop endif ; Fin de controle ; Suppression eventuelle des harm obtenus avec : ; . E0 (periode double en EW et distrib asym dans le plan (u,v)), ; . E1 avec les EW, ; . E2 (distrib asym dans le plan (u,v), ; . NS45 (barre en u pour v=45) ; . harmoniques obtenues avec des antennes trop proches de NS0. ; . suppression des doublures mediocres (bases NS courtes avec NS0 qui ; est parasitee, de toutes les bases EW avec NS0). ; . ; Voir description rangement dans harm_N dans IM_2D_CON. ; Rem (25 jan 02) : la mise a zero des harmoniques resoud le probleme ; de leur suppression : ; . s'il n'y a pas redondance, c'est evident. ; . s'il y a redondance il faut decider d'incrementer le compteur de ; redondance dans HARM_N_VIS (appele par IM_2D_CON, lui meme ; appele par CALC_IMAGE_HELIO) dans le cas de la nullite stricte) ; de la valeur des harmoniques. ; Proceder ainsi est plus facile que de definir des tableaux de per- ; mission d'utilisation d'antennes. ; Tests de la gestion des harmoniques. toto = moment(indgen(3)) ; pour compiler moment sans couper la liste. i_gestion = 0 if (i_gestion eq 1) then begin print, 'RH_IMAGE_CON , gestion des harmoniques :' if (correl(0, 18*25) eq 2) then begin print, ' Donnees : harm E0 * res NS, ' + $ 'pas d harm NS45 * res NS' i_e0_ns45 = 0 ; E0 correle avec les NS endif if (correl(0, 18*25) eq 43) then begin print, ' Donnees : harm NS45 * res NS, ' + $ 'pas d harm E0 * res NS' i_e0_ns45 = 45 ; E0 correle avec les NS endif if (ie0_ew eq 0) then print, ' Suppression harm E0 * res EW' if (ie1_ew eq 0) then print, ' Suppression harm E1 * res EW' if (ie0_ns eq 0) then begin ; E0 * res NS seulement avant jul_dat if (correl(0, 18*25) eq 2) then $ print, ' Suppression harm E0 * res NS' endif if (ie2_ew eq 0) then print, ' Suppression harm E2 * res EW' if (ie2_ns eq 0) then print, ' Suppression harm E2 * res NS' ; Tester si les harm de NS45 avec le reseau EW sont nuls ou non h45_ew = harm_N(432:449) mom = moment(h45_ew) c0 = complex(0, 0) if ( (mom(0) ne c0) and (mom(1) ne c0) ) then begin ch_45_1 = ' harm NS45 * res EW non nuls dans les donnees' ch_45_2 = ' ' endif if ( (mom(0) eq c0) and (mom(1) eq c0) ) then begin ; on admets que h45_ew est identiquement nul. Le if ne marche ; pas sur des vecteurs. ch_45_1 = ' harm NS45 * res EW nuls dans les donnees' ch_45_2 = ' (deja nuls dans les donnees)' endif if (ins45_ew eq 1) then print, ch_45_1 if (ins45_ew eq 0) then print, $ ' Suppression harm NS45 * res EW' + ch_45_2 ; Tester si les harm de NS45 avec le reseau NS sont nuls ou non h45_ns = harm_N(450:485) mom = moment(h45_ns) if ( (mom(0) ne c0) and (mom(1) ne c0) ) then begin if (i_e0_ns45 eq 0) then $ ch_45_1 = ' harm E0 * res NS non nuls dans les donnees' if (i_e0_ns45 eq 45) then $ ch_45_1 = ' harm NS45 * res NS non nuls dans les donnees' ch_45_2 = ' ' endif if ( (mom(0) eq c0) and (mom(1) eq c0) ) then begin if (i_e0_ns45 eq 0) then $ ch_45_1 = ' harm E0 * res NS ' if (i_e0_ns45 eq 45) then $ ch_45_1 = ' harm NS45 * res NS ' ch_45_2 = ' deja nuls dans les donnees' endif if (ins45_ns eq 0) then if (correl(0, 18*25) eq 43) then print, $ ' Suppression harm NS45 * res NS' if (ins45_ns eq 0) then print, $ ' Suppression harm NS45 * res NS' + ch_45_2 endif ; Fin des tests de gestion des harmoniques. ; Controle repris de RH_CALIB_CON icon = 0 if (icon eq 1) then begin print, 'RH_IMAGE_CON entree' hNS_NS0_ew = complexarr(25) for i=1, 24 do hNS_NS0_ew(i) = harm_N(9 + 18 *i) ch_format ="(50f8.3)" print, ' ' print, 'NS0_ew avec les NS' print, format=ch_format, abs(hNS_NS0_ew) print, ' ' print, 'NS0_ns avec NS7 a NS23' print, format=ch_format, abs(harm_N(573:575)), abs(harm_N(558:572)) print, ' ' print, 'NS8 avec NS0 a NS17' print, format=ch_format, abs(harm_N(540:557)) print, ' ' print, 'NS45 avec les NS' print, format=ch_format, abs(harm_N(450:485)) print, ' ' stop endif ; Fin de controle ; jul_ext = julday (10, 19, 1998) ; 19 oct 98. ext comme "extension". ; mise en service de NS45. if (ie0_ew eq 0) then harm_N (522:539) = 0 if (ie1_ew eq 0) then harm_N (504:521) = 0 if (ie0_ns eq 0) then begin ; E0 * res NS seulement avant jul_dat if (correl(0, 18*25) eq 2) then harm_N(450:485) = 0 ; la serie de NS45 avec les NS commence a 18*25 = 450. endif if (ie2_ew eq 0) then begin harm_N (486:503) = 0 ; E2 avec H1 a H16. harm_N (504) = 0 ; E1 avec E2. harm_N (522) = 0 ; E0 avec E2. endif if (ie2_ns eq 0) then for i=0, 24 do harm_N(18*i) = 0 if (ins45_ew eq 0) then harm_N(432:449) = 0 if (ins45_ns eq 0) then begin if (correl(0, 18*25) eq 43) then begin harm_N(450:485) = 0 ; l'ant 43 est NS45. Quand la NS45 est utilisee (a partir ; du 19 oct 98 environ), les harmoniques NS45*resNS sont ; a ces adresses, qui auparavant correspondaient a des cor- ; relations E0 * res NS. harm_N(572) = 0 endif endif ; Controle rtepris de RH_CALIB_CON icon = 0 if (icon eq 1) then begin print, 'RH_IMAGE_CON' hNS_NS0_ew = complexarr(25) for i=1, 24 do hNS_NS0_ew(i) = harm_N(9 + 18 *i) ch_format ="(50f8.3)" print, ' ' print, 'NS0_ew avec les NS' print, format=ch_format, abs(hNS_NS0_ew) print, ' ' print, 'NS0_ns avec NS7 a NS23' print, format=ch_format, abs(harm_N(573:575)), abs(harm_N(558:572)) print, ' ' print, 'NS8 avec NS0 a NS17' print, format=ch_format, abs(harm_N(540:557)) print, ' ' print, 'NS45 avec les NS' print, format=ch_format, abs(harm_N(450:485)) print, ' ' stop endif ; Fin de controle ; Suppression des petites bases a partir de NS0 (le "gap") if (gap_ns0(0) ne 0) then begin gew = gap_ns0(0) harm_N(9-gew:9+gew) = 0. endif if (gap_ns0(1) ne 0) then begin gns = gap_ns0(1) for j = 0, gns do begin harm_N(9 + 18*j) = 0. endfor endif ; Suppression des redondances if (i_redond eq 0) then begin ; usuel i_redond = 1 dans RH_CALIB_CON harm_N(6:14) = 0 ; bases EW a partir de NS0 de -400m a ; + 400m incluses. for j = 1, 6 do begin harm_N(9 + 18*j) = 0 ; NS0_ew avec NS1 a NS6 endfor harm_N(540:544) = 0 ; bases vers nord de NS8 avec NS0 a NS4 endif ; Suppression des harmoniques impliquant des antennes supposees a problemes harm_N_1 = harm_N(1) ; sauvegarde harmonique (NS0_ns, E1) ; pour le trac'e. d_0 = where(sel_ant eq 0, ic) ; numeros des antennes a eliminer. ; Reduction de sel_ant si les AA ne sont pas utilisees if (nbre_harm eq 576) then $ sel_ant = sel_ant(0:43) ; pour ne pas tenter d'annuler des harmoniques inexistants. if (ic ge 1) then begin print, 'RH_IMAGE_CON : mise a zero de tous les harm de ' + $ 'certaines antennes' for i=0, ic-1 do begin for j=0, nbre_harm - 1 do begin if ((correl(0, j) eq d_0(i)) or (correl(1, j) eq d_0(i))) $ then begin harm_N(j) = 0 ch_format = "(' harm_N(', i3, ') = 0')" ; print, format=ch_format, j endif endfor endfor endif ; Controle des harm NS avec NS0 (25 sep 03 icon = 0 if (icon eq 1) then begin harm_ns = complexarr(25) for j = 1, 24 do harm_ns(j) = harm_N(9 + 18*j) nz = 24 ; zoom de 0 a nz. NUMERO D'ANTENNE ET NON D'HARMONIQUE window, 0 plot, abs(harm_ns(0:nz)) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'abs(harm avec NS0_ew)' stop endif ; Controle de la distribution des harmoniques dans "harm_N". icon = 0 ; 1 pour mise au point. if (icon eq 1) then begin bmp = reform(harm_N, nbre_harm) help, harm_N, bmp mod_amp = abs (bmp) help, mod_amp plot, mod_amp(0:35) ; L'ordre de rangement des correlations est: correlations ; NS0 avec E2, E1, H1 a H16 0 a 17 ; NS1 ---- ---------------------------- 18 a 35 stop endif ; Fin de suppression eventuelle de certains harmoniques. ; Calcul de l'image (non polar ou polar) sur la grille des coord interferome- ; triques entieres. i_stop = 1 ; 1 pour arrets controles IM_2D_CON, $ ; clone de RH_MALC_IM_2D ; tous les parametres sont entree. sel_ant, $ i_stop, i_sys_lin, $ freq, ch_pol, ie0_ew, $ harm_N, harm_N_1, sref, correl, $ npi , n_bord, i_AA, $ icon_1, icon_2, icon_3, icon_4, $ n_win, frac_z, ch_siz, $ n_image, hmsc_deb, hmsc_fin end ; fin de RH_IMAGE_CON