PRO DP_CON_1, freq, ch_polar, harm_N ; Creation 10 avril 03 ; But : controler la correction des harmoniques pour la rotation et le depoin- ; tage fait pas CORR_ROT_DEP, a partir d'observation qui peuvent compor- ; ter des corrections fausses (CRE_CONFIG_1). ch_freq = strcompress(string(nint(freq))) + ' MHz ' h_NS0_ew = complexarr(24) h_NS0_ns = complexarr(24) for i=1, 23 do h_NS0_ew(i) = harm_N(9 + 18 *i) h_NS0_ns( 7) = harm_N(573 : 575) h_NS0_ns(10) = conj(harm_N(558 : 571)) ; pour avoir des bases orien- ; tees vers le nord dans les ; 2 cas. a_h_NS0_ew = abs(h_NS0_ew) a_h_NS0_ns = abs(h_NS0_ns) amax = max(a_h_NS0_ew) > max(a_h_NS0_ns) h_ew_1 = h_NS0_ew + 0.0001 * amax ; precaution pour atan. h_ns_1 = h_NS0_ns + 0.0001 * amax ; -------------------- p_h_NS0_ew = 180/!pi * atan(imaginary(h_ew_1), float(h_ew_1)) p_h_NS0_ns = 180/!pi * atan(imaginary(h_ns_1), float(h_ns_1)) ; Traces amplitudes amax = max(a_h_NS0_ew) > max(a_h_NS0_ns) amin = min(a_h_NS0_ew) < min(a_h_NS0_ns) da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da window, 0 plot, a_h_NS0_ew, xstyle=1, yrange=[y1, y2], ystyle=1 oplot, a_h_NS0_ns, linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_freq + ch_polar + $ 'abs(h_NS0_ew) ___, abs(h_NS0_ns) ....' ; Traces phases amax = max(p_h_NS0_ew) > max(p_h_NS0_ns) amin = min(p_h_NS0_ew) < min(p_h_NS0_ns) da = amax - amin y1 = amin - 0.05 * da & y2 = amax + 0.05 * da window, 1 plot, p_h_NS0_ew, xstyle=1, yrange=[y1, y2], ystyle=1 oplot, p_h_NS0_ns, linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, ch_freq + ch_polar + $ 'phase(h_NS0_ew) ___, phase(conj(h_NS0_ns)) ....' stop end ; fin de DP_CON_1 . ;______________________________________________________________________________ PRO RH_DPATCHFITS_NRH, $ ; entree. ; Caracterisation de l'observation nomfich, sel_freq, $ i_polar, t_deb, t_fin, $ ; Repertoire de sortie out_dir, $ ; Structure contenant les parametres de calcul d'image str_dp, $ ; sortie out_file ; Appel dans RH_DPNEW : ; RH_DPATCHFITS_NRH, $ ; entree. ; nomfich, sel_freq, $ ; i_polar, t_deb, t_fin, $ ; out_dir, $ ; str_dp, $ ; sortie ; out_file ;+ *********************************************************************** ; NAME: ; DPATCHFITS ; PURPOSE: ; Cette procedure cree des fichiers FITS 2D ; sur le repertoire out_dir ; -calcule les images heliographiques de np*np points sur un champ de ; largeur larg (RS). ; -selectionne les frequences a depatcher ; CATEGORY: ; Traitement de fichiers NRH ; CALLING SEQUENCE: ; RH_DPATCHFITS_NRH, $ ; nomfich, sel_freq, i_polar, t_deb, t_fin, out_dir, $ ; str_dp ; INPUTS: ; nomfich nom du fichier a traiter, ; i_polar 1 avec traitement des sorties polar, 0 sans. Si i_polar=1 pour ; une observation de calibration posterieure au 7 oct 02 les ; dipoles des antennes H2 H4 et NS0 doivent etre tournes de ; 45 degres. ; np dimension de l'image finale heliographique. ; mode_interpol interpolation mode from interferometric to heliographic image ; 1 linear, 2 cubic spline, 3 fft ; champ_helio largeur (4 ou 6 Rs) du champ heliographique (dans lequel le ; soleil est rond). ; t_deb heure de debut (h, m, s, c) ; t_fin -------- fin ----------- ; rot_p correction angle p : 0 sans, 1 avec (mais le centre du disque ; est toujours au centre du champ). ; out_dir directoire d'ecriture des fichiers FITS 2D ; sel_freq strarr(10) de selection des frequences ('Y' ou 'N') ; t_b choix des unites de brillance : 1 pour sfu, 2 pour Kelvins. ; COMMON BLOCKS: ; RH : communique avec les routines de lecture des donnees ; brutes RH ; EXAMPLE: ; ; MODIFICATION HISTORY: ; Ecrit par A. Bouteille et C. Mercier ; 97 sept - appel a INIT_MALAX ; - appels a maille_2d et ronron modifies ; - header FITS BINARY TABLE pour fichiers comprimes et decomprimes ; 97 oct - selection des frequences ; 98 aug 28 - appel des nouvelles routines read_rh (pour ajout antenne NS24) ; - attention a l'ordre des gains modifie dans readch_rh ; 99 jan 14 - correction des phases de Jean Marc ; avr 1 - MAILLE_2D est remplace par MAILLE_3. "maille" devient un ; fltarr(6) qui inclut les coord interfero du centre du disque ; (anciennement xc et yc). La maille calculee est la maille ele- ; mentaire. Donc ma/np est change en ma dans les appels ulteri- ; eurs. ; - RONRON, qui reprenait le calcul de la grille des coord interfe- ; ro pour calculer l'image ronde (ce que faisait aussi CALC_FLAG, ; remplace par GRILLE) est debarasse de ce calcul. Reduit a une ; ligne, RONRON disparait comme sous-programme. ; - DALC_IM_2D est epure et remplce par MALC_IM_2D ; 99 mai 6 - reconstitution approximative de la composante continue, pour ; l'image et le lobe theorique, comme moyenne des bas harmoni- ; ques. Permet de normaliser l'image a celle obtenue avec un lobe ; d'integrale 1 en utilisant directement le lobe theorique (aupa- ; ravant depourvu de composante continue). ; 99 aug 22 - le calcul de "champ", drapeau de limitation du calcul du champ, ; est corrige et reporte dans GRILLE. Il tient compte: ; . du champ max demande par l'utilisateur, ; . du champ interferometrique qu'il ne depasse pas. ; 99 aug 26 - appel separ'e de INIT_MALAX_2D (nouvelle version d'INIT_MALAX ; qui prend ie0 en entree au lieu de le mettre toujours a 1) par ; INIT_MALAX_RH_2D (nouvelle version de INIT_MALAX_RH qui passe ; a INIT_MALAX_2D la valeur de ie0 choisie dans DPATCHFITS). ; 99 sep 13 - drapeaux d'utilisation en entree de DPATCHFITS: ; . ie0_ns pour harm E0 avec les NS ; . ie0_ew --------------------- EW ; . ie2_ns --------- E2 -------- NS ; . ie2_ew --------------------- EW ; 99 oct 7 - modularisation, regroupement dans calc_image_helio ; 8 - calibration par methode des antennes relais. ; 99 oct 18 - definition d'une zone d'exclusion autour de NS0 (harmoniques ; trop parasites). ; 00 dec 1 - nouveaux arguments en entree : ; mode_clean 0 pour clean normal sur champ entier, ; 1 sur champ limit'e (suppression de centre d'orage) ; x_dom, y_dom, coord (RS, >0 ou <0) centre du domaine de clean. ; dr_dom, largeur totale (RS) du doamine de clean. ; dec 14 - n_period est defini dans dpnew et passe en argument a DPATCHFITS ; 01 fev 12 - 4 modes de clean (voir mode_clean dans dpnew). ; . 0 pas de clean, ; . 1 pour clean normal sur champ entier, ; . 2 pour chercher partout les comp sales mais ne pas restituer ; les composantes synthetiques interieures a dom_c_i ; defini par x_dom, y_dom et dr_dom. ; . 3 pour retirer les composantes sales de dom_c_i et obtenir ; une image residuelle non nettoyee du reste. ; . 4 pour nettoyer (rechercher les composamtes sales et resti- ; tuer les composantes synthetiques dans dom_c_i seulement. ; - definitions de tous les parametres de description du clean dans ; dpnew. ; 01 jun 8 - introduction de i_sys_lin pour choisir la methode de resolution ; des systemes lineaires : 1 Cramer, 2 LU, 3 SVD. Placee en ; amont de CALIBRATION.PRO et SIMUL_HARM_N, la modifica- ; tion de i_sys_lin se fait en un seul endroit. ; 01 jul 6 - definition dans DPATCHFITS du nbre de points npi=128 de l'image ; interferometrique. On ne prend plus npi = np . ; sep 10 - definition de ie0, ie1, ie2, ins45 transferee de DPATCHFITS dans ; INIT_MALAX_2D. ; 19 - le choix du mode d'interpolation pour passer de l'image interfe- ; rometrique a l'image heliographique est remont'e dans DPNEW, ; compte tenu que les fft sont lentes. ; 20 - on remonte par CALC_IMAGE_HELIO le flux total et le flux compact ; calcules par MALC_IM_2D. ; 25 - variable i_pave_plein passee a SIMUL_HARM_N et CALIBRATION, ; CALC_IMAGE_HELIO pour MALC_IM_2D. ; oct 18 - par suite de l'abandon de toute forme de delouiserie il devient ; inutilde de faire figurer i_filtr_ech (l'ancien i_choix) dans ; la liste d'appel. La valeur de i_ech_max donne l'information ; equivalente. ; 02 fev 16 - ajout de gain_ant et phase_ant) dans la liste de sortie de ; SIMUL_HARM_N pour passage a CALIBRATION pour comparaison ; avec amul et bmul lors de la verification finale. ; 02 mar 5 - ajout de sel_ant dans la structure des parametres d'entree pour ; eliminer les antennes a problemes. ; mar 7 - ajout de crit_arret. larg_2 devient un intarr(2) ; 19 - reorganisation : la calibration directe ou auto) est faite apres ; integration des donnees : on sort la calibration de la boucle ; de traitement des scrutations. ; - la calibration peut etre directe (amorcage avec un modele de ; Cygne simple). ; 27 - introduction de rot_p pour corriger ou non de la rotation d'angle ; p. ; 28 - passage de t_moy a SIMUL_HARM_N poir simuler un Cygne reel a une ; heure choisie (la maille de la grille de harm de harm_N n'est ; alors plus carree). ; 02 oct 28 - extraction de la declinaison, du depointage du reseau H et de la ; rotation de H2, H4 et NS0 pour transmission s CALIBRATION pour ; pouvoir utiliser les anciennes observation "sources de contro- ; le" pour recalibrer a posteriori. ; dec 5 - mise de calib_np_p dans str_dp ; - introduction de "out_file" en sortie de RH_DPATCHFITS_NRH ; 03 mar 13 - ajout de i_recentrage a la liste d'entree pour passage a ; RH_MALC_IM_2D a travers RH_CALC_IMAGE_HELIO ; 03 mar 18 - elimination de la procedure de calibration, deja regoupee dans ; RH_DP_CALIB ; 03 mar 19 - ajout de x_centrage et y_centrage (en canaux interferometriques) ; a la liste d'entree pour mieux desaliaser. ; 03 mar 20 - tf_objet_simule transmis de SIMUL_HARM_N a RH_CALC_IMAGE_HELIO ; 03 jun 23 - ajout des gains et phases de E2, E1, E0, NS0_ew, H15, H16, ; NS0_ns, NS8m NS21, NS22, NS23, NS45 en entree pour correction ; pragmatique de l'instabilite du reseau. ; 03 jun 24 - les parametres pour la recalibration et/ou l'auto-re-calibration ; partielle, ainsi que les parametres de recentrage pour remplir ; le pave central sont incorpores dans la structure str_dp. ; 03 sep 17 : ajout a la structure str_dp : ; . i_corr_pos_ant, i_source_cal, h_cal, ; retrait de calib_np_p ; 04 jan 12 : ajout de g_AA1 etc aux liste des *STR_DP ; 04 sep 22 : refonte du systeme de correction manuel des antennes : g_E2 etc. ; disparaissent et on introduit les tableaux(48) g_flag, g_amp, ; g_phi. ; 05 feb 2 : ajout de fac_integ, nbre de scrutations cumulees avant traitem- ; ent, a la liste d'entree de PUT_STR_DP, GET_STR_DP, DEF_STR_DP, ; HDR_STR_DP. Cela permet d'integrer souplement un fichier d'en- ; tree a haute cadence, par exemple pour avoir une cadence de ; sortie comparable a celle du GMRT. ; 05 jun 15 : x_centrage, y_centrage (definis en 0.01 Rs dans RH_DP_NEW) sont ; convertis en canaux interferometriques dans RH_DPATCHFITS_NRH ; au debut du "traitement de calcul d'image" (on ne les utilise ; que dans MALC_IM_2D, mais il faut disposer de l'heure hmsc). ; 05 nov 24 : chgt du while pour pouvoit traiter une seule scrutation, p ex ; dans des fichiers inegres 10 sec ou 128 sec : ; while (klu le kfin) -> while (klu le kfin + 6) ; l'ajout de 6 est destine a taiter le cas ou on veut isoler une ; seule scrutation (dans un fichier 10sec ou 128 sec), et entrer ; dans la boucle si une une frequence autre que la frequence 0 ; est selectionnee. Dans ce cas on peut prendre t_fin = t_deb. ;************************************************************************** ; NOTATIONS: ; entfi.freq frequences (intarr(10)) ; entfi.corel tableau (2, 576 ou 648) indiquant les numeros des 2 antennes ; correspondant a chacun des 576 (ou 648) harmoniques. Permet ; dans MALC_IM_2D le passage de harm_N au tableau des harmoni- ; ques a 2 dimensions. ; bfi.h (h,m,s,c) de debut d'observation. ; klu, kdeb, kfin, kmax numeros des enregistrements. ; numim numero d'image. ; sfu_tb ; lufi.pt fltarr(4, 576 ou 648) ; i_redond 0 pour supprimer les doublures les moins bonnes. ; gap_ns0 intarr(2), definit une zone d'exclusion d'harmoniques parasi- ; te's autour de NS0. ; sel_ant drapeau d'utilisation des antennes (0 si antenne a probleme). ; entfi.frq tableau des frequences de l'observation ; g_NS45 tableau : drapeau, gain et phase (deg) de NS45. ; g_AA1 tableau : drapeau, gain et phase (deg) de AA1. ; isel tableau des frequences qu'on a choisi de traiter ; nsel nombre de frequences a traiter ; entfi.nval integer (2304 par exemple). ; pt fltarr(4, 576 ou 648) puis fltarr(2, 2, 576 ou 648). Tableau ; des harmoniques : ; 1er indice: cos ou sin, ; 2eme ind: non polar ou polar, ; 3eme indice: balaie les 576 ou 648 harm. ; xn, yn fltarr(np, np). Coordonnees heliographiques normalisees des ; points de l'image ou le soleil est rond. Entre -1 et 1-1/64. ; champ_helio largeur totale (Rs) du champ couvert par la grille heliographi- ; que a mailles carrees (dans laquelle le soleil est rond). A ; ne pas confondre avec la largeur du champ interferometrique. ; np dimension demandee pour l'image heliographique finale (soleil ; rond dans un champ de largeur "larg" RS). ; mode_clean choix du mode de clean : ; . 0 pas de clean ; . 1 pour clean normal sur champ entier, ; . 2 pour chercher partout les comp clean mais ne pas resti- ; tuer celles interieures au domaine dom_c_i defini par ; x_dom, y_dom et dr_dom. ; . 3 pour retirer les composantes clean de dom_c_i et obtenir ; une image residuelle non nettoyee du reste. ; . 4 pour nettoyer seulement dans dom_c_i . ; x_dom, y_dom coord (RS, >0 ou <0) centre du domaine de clean. ; dr_dom, rayon (RS) du domaine de clean. ; dom_c_i tableau des indices de l'image interferometrique du domaine de ; clean ; DEBUT EFFECTIF DE LA PROCEDURE version_solarsoft = 1 ; 0 version test de Claude 1 SolarSoft ; permet de gerer impressions et stop RH_CORR_ANTENNES ; Compilation de : ; CORR_POS_ANT correction de position des antennes ; CORR_GAIN_ANT correction manuelle du gain de certaines antennes ; CORR_ROT_DEP_NP correction depointage et rotation non polar ; CORR_ROT_DEP_P --------------------------------- polar if version_solarsoft eq 0 then $ RH_GENE_SUB ; Compilation de : ; DISQUE ; ANNEAU (et d'autres, inutiles ici) ; Restauration des variables contenues dans la structure str_dp GET_STR_DP, str_dp, $ version, i_simulation, $ ; Temps d'integration prealable t_integ, $ ; Caracterisation de l'image n_period, mode_interpol, rot_p, np, champ_helio, t_b, $ i_affich_lobe, i_mul_pol, $ ; 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, $ ; Cas particulier des sources de controle de calibration i_con_cal, $ ; Correction de position d'antennes et recalibration i_corr_pos_ant, i_recalib, i_source_cal, h_cal, ch_cor, $ phi_ew_ns, phi_w_e, i_corr_g, g_flag, g_amp, g_phi, $ ; Complement du pave central (harm multiples impairs de 50m) i_pave_plein, i_recentrage, x_centrage, y_centrage, $ ; Calibration et clean mode_clean, x_dom, y_dom, dr_dom, $ i_ech_min, i_ech_max , a_in, a_ex, $ a_red_fil, itermax, gain_iter, jmax, crit_arret, $ i_seuil, i_test_arret, i_stop_image, i_ecrete, frac_ecr, $ tete_fich, $ fac_crois_res, coeff_seuil, frac_flux, $ i_lissage, poids_c, larg_2 x_centrage_0 = x_centrage y_centrage_0 = y_centrage @rh_common.inc ; inclut le contenu de "common_rh.inc" (description de ; la structure "entfi." (en-tete de fichier). ; Voir dans : ; /home/dirhelio/sswnrh/radio/nrh/idl/-> ; dp_soft_rh/routine_rh ; et faire xemacs rh_open.pro ; STRUCTURE entfi (entete), repris de de rh_common.inc ; ======================== ; ; Longueur tete du fichier entFI.lgent Long ; Longueur du fichier entFI.filesize Long ; Longueur d'une acquisition, en octets entFI.lg Int ; Type des acq du fichier, 1-1D, 2-2D entFI.typ Int ; Codage des acquisitions, 2-16 bits entFI.cod Int ; Representation des acq, 0-Harm entFI.rep Int ; Nombre de valeurs tete acq entFI.lgtet Int ; Nombre de valeurs corps acq entFI.nval Int ; Nombre de valeurs 2D entFI.nval_2d Int ; Nombre de valeurs EW entFI.nval_ew Int ; Nombre de valeurs NS entFI.nval_ns Int ; Type de donnees,'1D ','2D ' entFI.d_typ Bytarr(6) ; Temps d'integration en millisecondes entFI.itg Long ; Date (jour, mois, an) entFI.dat Intarr(3) ; Reference de la source, en caracteres entFI.ref Bytarr(6) ; Heure meridien (Heure, minute, sec, centieme) entFI.mer Intarr(3) ; Declinaison (degres,',',') entFI.dec Intarr(3) ; Horloge 0-TU ou 1-TS entFI.hg Int ; Heure de debut (heure, minute) entFI.hdeb Intarr(4) ; Heure de fin en heure, minute) entFI.hfin Intarr(4) ; Nombre de frequences, de 1 a 10 entFI.nf Int ; Table des frequences en centaines de Khz entFI.frq Intarr(10) ; Correction de trajectoire 1-avec ou 0-sans entFI.trj Int ; Compression 1-avec ou 0-sans entFI.comp Int ; Cycle recepteur, en millisecondes entFI.cyclms Int ; Description donnees entFI.d_obs bytarr(78) ; Reseau EW depointage 0-Sans, 1-Est, 2-Ouest entFI.depew Int ; Reseau EW depointage en minutes entFI.dephew Int ; Rotation dipoles ant 2, 4 EW 1-avec ou 0-sans entFI.rotew Int ; Rotation antenne 0 NS 1-avec ou 0-sans entFI.rotns Int ; Position EW des antennes en mm entFI.pxant Lonarr(44) ; Position NS des antennes en mm entFI.pyant Lonarr(44) ; Table des correlations (no Ref, no ant) entFI.corel intarr(2,576) ; Tables de correction entCOR:entCOR ; Description de l'observation descrip bytarr(80,2) ; Description de l'activite activ bytarr(80) ; Description des pannes pannes bytarr(80,3) ; Evenements non comprimes evenements bytarr(80,15) ; swap_endian a faire si = 1 endian 0L ; numero du dernier enregistrement klumax 0L fich_cor = ch_cor ; en attendant les chgts de notation dans les program- ; mes *str_dp (c'est juste une question de notations, ; ca transite par str_dp). ; Annulation d'ensemble de la correction de certaines antennes g_flag = i_corr_g * g_flag ; Verification de non nullite de ces gains for ia = 0, 47 do begin if (g_amp(ia) eq 0) then begin print, 'Erreur : un des gains des antennes est nul.' print, 'Taper .c pour continuer quand meme.' stop endif endfor ; Choix de la methode de resolution des systemes lineaires apparaissant dans ; la determination de parametres aux moindres carres dans fit_quadric, ; fit_cubic, fit_quartic (eux memes appelles par malc_im_1d et ; malc_im_2d (contenus dans CALIB_SUB) et SOUSTRAC_BASE. i_sys_lin = 2 ; 1 Cramer, 2 LU, 3 SVD. npi = 128 ; nbre de points sur l'image interferometrique (ce n'est plus ; np). ; Ouverture du fichier de type Nancay ; if (not open_rh(nomfich, /mono, /malax)) then begin ; jusqu'au 26 aug 99. if (not RH_OPEN(nomfich, /mono )) then begin ; apres le 26 aug 99. print, 'fichier de donnees non trouve',nomfich return endif ; Le if (...) appelle open_rh (avec le seul mot-cle mono) puis teste ; les codes de sortie. S'ils ne sont pas bons, impression d'un code ; erreur. ; Tests de securite (a completer) (a mettre apres le RH_OPEN, pour avoir la ; structure "entfi" de l'observation demandee, et non celle correspon- ; dant a un precedent passage) if ( (t_deb(0) eq t_fin(0)) and $ (t_deb(1) eq t_fin(1)) and $ (t_deb(2) eq t_fin(2)) and $ (t_deb(3) eq t_fin(3)) ) then begin print, '' print, 't_fin = t_deb => on prends t_integ = -1 (pas d''inte' + $ 'gration) et on continue' ; stop endif if ((t_integ ne -1) and (t_integ le 0.001 * entfi.itg)) then begin ; car t_integ est en secondes et entfi.itg est en msec. print, '' print, 'Le temps d''integration demand''e t_integ est trop ' + $ 'court pour le fichier de' print, ' donnees' print, '' print, 'Taper .c pour prendre t_integ = -1 (pas d''integration) ' + $ 'et pour continuer' print, '' stop t_integ = -1 print, format="('nouvelle valeur de t_integ =', i3)", t_integ print, '' endif ; La structure str_dp_hdr contient les parametres du calcul a afficher ; dans le header fits. C'est un sous ensemble de str_dp. HDR_STR_DP, str_dp, $ ; entree str_dp_hdr ; sortie ; Calcul du facteur d'integration if (entfi.comp eq 0) then begin ; donnees non comprimees dom = where(sel_freq eq 'Y', ic) if (ic eq 1) then begin ; une seule freq selectionnee. if (t_integ le 0) then fac_integ = 1 if (t_integ gt 0) then fac_integ = $ nint(t_integ / (0.001 * entfi.itg)) endif if (ic ne 1) then begin ; plusieurs freq selectionnees fac_integ = 1 print, 'Pas d''integration pour plus d''une frequence ' + $ 'selectionnee' if version_solarsoft eq 0 then begin print, 'Taper .c pour continuer, sans integration' stop endif endif endif ; fin du if donnees sans compression. if (entfi.comp ne 0) then begin ; donnees comprimees fac_integ = 1 endif correl = entfi.corel ; alleger notations et retablir orthographe. taille = size(correl) nbre_harm = taille(2) ; Lecture de l'heure de calibration dans l'entete et ecrasement de h_cal ; On prend en fait l'heure de calibration correspondant a la 1ere ; frequence selectionnee, admettant que c'est la meme pour toutes les ; frequences. dom = where(sel_freq eq 'Y', ic) if (ic eq 0) then begin print, 'Aucune frequence selectionnee. STOP !' endif h_cal_en_tete = entfi.entcor.heure ; Si on fait directement : ; h_cal_en_tete = reform(entfi.entcor.heure(*, dom(0))) ; ca laisse un intarr(2, 10). On commence donc par reduire au statut ; de variable normale h_cal_en_tete = reform(h_cal_en_tete(*, dom(0))) ; Ecrasement de h_cal pour les observations comportant l'eure de calibration if ((h_cal_en_tete(0) ne 0) and (h_cal_en_tete(1) ne 0)) then $ h_cal = h_cal_en_tete ; Initialisation du common MALAX ; Ce common est utilise par : ; . RH_GRILLE, RH_MAILLE_3, CALCPOS, pour les soleil, ; . VISIB_CYGNE pour une calibration directe. ; Drapeaux d'utilisation des antennes d'extension EW ; L'acquisition des donnees a Nancay s'est faite de deux facons: ; - jusqu'au 7 oct 98 inclus l'antenne de reference pour les correla- ; tions 453 a 485 est E0 (avec les NS), ; - a partir du 8 oct 98 inclus E0 est remplacee par NS45 dans ces ; memes correlations. ; ; Dans la version de synthese des images 2-dim de Kerdraon seuls les ; harmoniques du "pave central" sont utilises. Les correlations de ; E0 et E2 et le champ interferometrique en EW (avec pas EW de 100m) ; etait insuffisant au dessus de 327 MHz (aliasing des bords E et W ; du soleil. On calculait deliberement l'image en coordonnees inter- ; ferometriques entieres sur un champ double du champ interferometri- ; que EW pour bien visualiser les effets de l'aliasing. ; Dans la version de synthese des images 2-dim de J-M toutes les cor- ; relations disponibles sont utilisees. Du fait de E0 le champ EW ; est, stricto sensu, double du precedent, mais ces harmoniques sup- ; plementaires etant tres minoritaires, cela ne remedie pas vraiment ; a l'aliasing. Cependant il suffit maintenant d'un champ egal au ; champ EW theorique pour le visualiser. ; Jusqu'au 7 oct 98 les harm impliquant E0 et E2 avec les NS ont la ; disposition: ; E2 E0 ^ v ; * * | ; * * | ; *-------*-------o-------*-------*-> u ; * * ; * * ; E0 E2 ; Apres le 7 oct 98 la "barre" decentree de E0 est remplacee par une ; "barre" de NS45 sur l'axe v. Ainsi les harmoniques avec une base EW ; multiple impaire de 50m sont encore plus minoritaires et l'alia- ; sing est encore moins reduit. Par contre la couverture (u, v) est ; meilleure. ; Apres le 16 sept 99 on peut choisir ie0_ew, ie0_ns, ie2_ew, ie2_ns. ; Initialisation du common MALAX, qui servira a passer des ccord interfero aux ; coord helio. POS_HELIO (appele par RH_CALCPOS) a besoin de connaitre ; le champ interferometrique et le nbre de pts sur ce champ. Ici on ne ; fixe pas le champ interferometrique, on ne fait qu'etre coherent ; avec le choix fait plus bas dans MALC_IM_2D, ou la place est prevue ; dans le tableau des harm a 2 dim pour les harm avec E0, meme s'ils ne ; sont pas utilises. Ainsi le champ interferometrique obtenu par fft ; correspond a un pas EW de 50m (ce qui fait 2 periodes EW si ; ie0_ew = ie0_ns = 0). ; Definition de variables necessaires pour calculer le common MALAX hmer = fltarr(3) hmer(0:2) = float (entfi.mer(0:2)) hmer(2) = hmer(2) + float(entfi.mer(3)) / 100. dec = fltarr (3) dec(0:2) = float (entfi.dec(0:2)) dec(2) = dec (2) + float(entfi.dec(3)) / 100. ; Declinaison, depointage, rotation des antennes, dephasage polar ; (utile pour la calibration) delta = !pi/180 * (dec(0) + float(dec(1)) / 60 + float(dec(2)) / 3600) ; les 0.01" d'arc sont inclus dans dec(2) rotew = entfi.rotew rotns = entfi.rotns if (entfi.depew eq 0) then ah_dep = 0 if (entfi.depew eq 1) then ah_dep = - entfi.dephew ; en minutes de temps. if (entfi.depew eq 2) then ah_dep = entfi.dephew ; ------------------- ah_dep = !pi / 12 * ah_dep / 60 ; en radians, >0 pour depointage W. date_obs = entfi.dat ; pour eviter date(0) qui plante (nouvelle ; version IDL ?). sref = string (entfi.ref) ; 'soleil' ou 'autre'. if (i_simulation gt 0) then begin ; On cherche a dependre le moins possible du fichier d'observation. hmer = [11, 50, 00] dec = [20, 0, 0] t_deb = [11, 50, 00, 00] ; usuel : 11, 50, 00, 00 t_fin = [11, 55, 00, 00] ; ----- 11, 55, 00, 00 endif ; Calcul du temps t_moy de milieu d'observation ; Dans une calibration directe, t_moy est utilise par ; . VISIB_CYG (dans CALIB_SUB) pour le calcul de la visibilite du mode- ; le de Cygne rustique. Le meme calcul aura eventuellement au debut ; pour simuler l'observation. ; . COMPAR_ORIG_RECAL_1D et COMPAR_ORIG_RECAL_2D (dans CALIB_SUB) ou ; il figure dans la liste d'entree mais est en fait inutilise. ah_deb = t_deb(0) + t_deb(1)/60. + t_deb(2)/3600. + t_deb(3)/360000. ah_fin = t_fin(0) + t_fin(1)/60. + t_fin(2)/3600. + t_fin(3)/360000. t_moy = (ah_deb + ah_fin) / 2 ; Initialisation du common MALAX INIT_MALAX_2D, date_obs, hmer, dec, sref, rot_p ; entree. ; - Appelle IN_MALAX_2D: initialisation du common MALAX, utilis'e ; par RH_CALCPOS (appele par RH_MAILLE_3 et ...) et qui contient: ; . date, nature source (soleil ou non), drapeaux de corr iono, ; . parametres des reseaux, ; . new et nns, nbre de canaux "minima" en EW et NS, ; . heure meridien (TU ou TS decimale), declinaison, ; . angle P et derive en H et delta (nuls pour sources non solaires) ; Notations ; entfi.dat date (j, m, a) ; hmer temps meridien (TU pour soleil, TS pour sources), ; c'est un fltarr(3), comme dec ; sref 'soleil' ou 'autre' ; rot_p correction de l'angel p : 0 sans, 1 avec. ; - la sortie est le common MALAX, utilis'e dans ; . RH_MALC_IM_2D pour convertir x_centrage et y_centrage ; (definis en 0.01 Rs) en canaux interferometriques, ; . RH_GRILLE, ; . RH_MAILLE ; ; common MALAX, $ ; ij, im, ian, icorpoi, icorion, isour, $ ; etmer, edec, dew, new, hew, u, dns, hns, a1, nns, sinl, cosl, $ ; sinp, cosp, gdel, ghmer, ahmer, c, rsol, ahmersol, sind, cosd ; Initialisation d'une correction de phase a faire jusqu'au 5/11/1998 non in- ; clus (correction d'un sabotage delouisien traitre, sorte de chauuse- ; trappe non signalee). ; treux). jul_cor = julday (11, 5, 1998) ; pour "julien correction". 5 nov 98 jul_dat = julday (date_obs(1), date_obs(0), rh_date100(date_obs(2), /full)) ; pour "julien date". icorrec_phase = 0 if (jul_dat lt jul_cor) then begin icorrec_phase = 1 print, 'Correction des phases: entfi.nval =', entfi.nval / 4 endif ; Conversion des temps de debut et de fin en msec temps = [3600000L, 60000L, 1000L, 10L] msd = (transpose(temps) # long(t_deb)) (0) ; h debut en msec msf = (transpose(temps) # long(t_fin)) (0) ; - fin ------- if (i_simulation gt 0) then begin nbim = 1 kdeb = 0L kfin = 0L klu = 0L np = 128 isel = where (sel_freq eq 'Y') nsel = n_elements (isel) numrec = lonarr(entfi.nf) numrec(*) = 0L numcnt = lonarr(entfi.nf) numcnt(*) = 0L iunit = lonarr(entfi.nf) imsdeb = 0L imsfin = 24L * 3600000L imstd = lonarr (entfi.nf) imstf = lonarr (entfi.nf) imstd(*) = imsdeb imstf(*) = imsfin itd = intarr (4, entfi.nf) itf = intarr (4, entfi.nf) mmax = fltarr (entfi.nf) mmin = fltarr (entfi.nf) hmsc = t_deb numtot = 1 goto, a0 endif ; Liste de variables : ?? (27 nov 01) ; Selection des heures: ; Calcul des indices des records de debut et de fin nbim = -1 if (msd eq msf) then nbim = 1 ; Donnees non comprimees if (entfi.comp eq 0) then begin ind = where (t_deb ne t_fin, count) ; nbre d'elements differents ; dans t_deb et t_fin. RH_HDEB, t_deb, $ ; entree, heure de debut demandee. kdeb, h_lu ; sortie, kdeb numero d'image de debut, suivant ; immediatement le debut demand'e ; t_deb, ; h_lu heure --------------- t_deb = h_lu ; heure de debut effective. print, 'numero enregistrement de debut=', kdeb ch_format = "('Heure trouvee juste apres le debut demande : ', " + $ "i2.2, ':', i2.2, ':', i2.2, ':', i2.2)" print, format= ch_format, h_lu if (count ne 0) then begin ; les heures de debut et fin demandees ; sont differentes. print, 'Heures de debut et de fin demandees differentes' RH_HDEB, t_fin , $ ; entree, heure de fin demandee. kfin, h_lu ; sortie, kfin numero d'image de fin, suivant ; immediatement la fin demandee t_fin ; h heure lue suivant immediatement ; l'heure de fin demandee. ch_format = "('Heure trouvee juste apres la fin demandee :" + $ " ', i2.2, ':', i2.2, ':', i2.2, ':', i2.2)" print, format= ch_format, h_lu print, ' (adoptee comme fin effective) t_fin = h_lu ; heure de fin effective. end else begin print, 'Heures de debut et de fin demandees identiques' kfin = kdeb t_fin = t_deb print, 'Heure de fin effective adoptee identique a l''heure' + $ ' de debut effectif' endelse end else begin ; Donnees comprimees RH_LECDEB, msd, msf, $ ; entree (heure debut et fin en msec). kdeb, kfin ; sortie numeros extremes des enrg a traiter. if (nbim eq 1) then kdeb = kfin endelse ch_format = "(8x, i6, 3x, i6, 6x, i6)" print, format=ch_format, kdeb, kfin, entfi.klumax ; Si klumax n'est pas une frequence selectionnee, on demande a lire une frequen ; ce au dela de klumax et on tombe sur "end of file". isel = where (sel_freq eq 'Y') nsel = n_elements (isel) ; if ( abs(kfin - entfi.klumax) lt entfi.nf) then $ ; kfin = entfi.klumax -entfi.nf + isel(nsel-1) + 1L print, 'RH_DPATCHFITS_NRH, numeros d''enregistrements (frequence ' + $ 'numero 0)' print, ' pour les heures de debut et de fin demandees :' print, ' debut fin fin fichier' ch_format = "(8x, i6, 3x, i6, 6x, i6)" print, format=ch_format, kdeb, kfin, entfi.klumax klu = kdeb if ( nsel eq 1) then begin ; Calcul des numeros d'enregistrement de debut et de fin pour la frequence ; selectionnee (et non plus la frequence 0) klu_0 = klu ; pour la frequence numero 0. RH_READCH, klu, sel_freq ; En entree sel_freq indique les frequences selectionnees. ; En sortie les variables sont dans le common RH. ; RH_READCH boucle sur les scrutations en incrementant klu, jusqu'a ; trouver une scrutation pour une frequence selectionnee. A ce moment ; on sort sans incrementer klu, qui doit etre increment par l'utili- ; sateur. ; Au final klu ainsi incerment'e en sortie de RH_READCH est le numero ; de scrutation pour la frequence selectionnee a l'heure ; qui suit immediatement l'heure de debut demandee. ; klu est ainsi aygment'e de 1, 2, 3 (par rapport a sa valeur initi- ; ale pour la frequence 0) pour les frequence de numeros 1, 2, 3. dklu = klu - klu_0 ; incrementation de klu pour trouver la ; frequence demandee. kdeb = kdeb + dklu ; pour la frequence selectionnee. ; pb si fin de fichier kfin = kfin + dklu ; ------------------------------ endif ; Si klumax n'est pas une frequence selectionnee, on demande a lire une frequen ; ce au dela de klumax et on tombe sur "end of file". if ( abs(kfin - entfi.klumax) lt entfi.nf) then $ kfin = entfi.klumax -entfi.nf + isel(nsel-1) + 1L if version_solarsoft eq 0 then print, ' ' ch_freq = string(format="(i3)", entfi.frq(isel(0))/10) + ' MHz' ; rappel : isel est le tableau des numeros de frequences selectionnees. ; isel(0) et la plus basse. print, 'Frequence ' + ch_freq + ' :' print, ' numeros d''enregistrements pour les heures de debut et de' + $ 'fin demandees :' print, ' debut fin' ch_format = "(16x, i6, 3x, i6, 6x, i6)" print, format=ch_format, kdeb, kfin klu = kdeb RH_READCH, klu, sel_freq ; Lecture prealable, faite ici seulement pour avoir l'heure de debut. hmsc = bfi.h ch_format = "('RH_DPATCHFITS_NRH, heure debut trouvee: ', " + $ "i2.2, ':', i2.2, ':', i2.2, ':', i2.2)" print, format=ch_format, hmsc ; Mise en forme des header fits, ouverture des fichiers et ecriture des ; headers provisoires. numrec = lonarr(entfi.nf) numrec(*) = 0L numcnt = lonarr(entfi.nf) numcnt(*) = 0L iunit = lonarr(entfi.nf) imsdeb = 0L imsfin = 24L * 3600000L imstd = lonarr (entfi.nf) imstf = lonarr (entfi.nf) imstd(*) = imsdeb imstf(*) = imsfin itd = intarr (4, entfi.nf) itf = intarr (4, entfi.nf) mmax = fltarr (entfi.nf) mmin = fltarr (entfi.nf) dirdsmf = strarr(entFI.nf) itg_lu = entfi.itg RH_INIT_HEADER, $ ; Entree $ nomfich, entfi.nf, entfi.frq, entfi.dat, bfi.h, $ kdeb, kfin, itg_lu, fac_integ, sel_freq, $ i_polar, np, $ imstd, imstf, champ_helio, t_b - 1, mmax, mmin, $ out_dir, str_dp_hdr, $ ; Sortie dsmf, numcnt, iunit, dirdsmf, errmsg ; Rem : on passe t_b - 1 a partir du 6 jul 01 (au lieu de ; l'ancien tb qui valait 0 ou 1 pour eviter de modifier ; D2D_FITS. ; Notations : ; dsmf nom de fichier de sortie ; numcnt tableau (dimensionn'e au nbre de frequences), qui ; donne le nbre d'images a ecrire pour chaque fre- ; quence. Les images en sortie sont comprimees si ; les donnees en entree le sont. if errmsg ne '' then return ; print,' sortie de RH_INIT_HEADER' out_file=dirdsmf ; Calcul du nombre total d'images a ecrire (toutes frequences jointes) numtot = 0L for i=0, entfi.nf-1 do $ ; somme sur les frequences traitees. if (sel_freq(i) eq 'Y') then numtot = numtot + numcnt(i) a0: J_M_Delouis = 0 ; Borne de fin du saut dans le cas de donnees simul'ees ; commemorant aussi la presence a cette place d'un ; fragment de calcul delouisien (calcul de xn et yn). numim = lonarr(entfi.nf) numim(*) = 0L ; Compilation preliminaire de RH_SIMULATION (aucune execution) qui contient ; SIMUL_HARM_N et SIMUL_RH (appele par MALC_IM_2D pour produire des ; images theoriques du soleil avec les antennes anti-aliasing. if (i_simulation gt 0) then RH_SIMULATION ; Boucle sur les acquisitions (1 acq = 576 ou 648 harmoniques) ; Les enregistrement des differentes frequences a une meme heure sont ; contigus (numero klu). La boucle sur les frequences est donc interne ; a la boucle sur les heures. ; Rappel : isel tableau des frequences qu'on a choisi de traiter. ; nsel nombre de frequences a traiter. ; nf numero de frequence traitee dans la boucle du while. n_image = 1 ; initialisation du numero d'image traitee. n_bord = 10 ; hauteur des bandes N et S pour ajuster base dans ; RH_MALC_IM_2D. Usuel 10. C'est suffisamment faible ; pour que le soleil n'y ait pas de contribution ; importante. Pour les sources non solaires c'est un ; peu faible pour determiner une bonne base, mais ; la limitation du domaine d'integration du flux ; dans RH_MALC_IM_2D (modif du 2 dec 2005) limite ; bien les erreurs sur le flux. ; print, 'Entree dans la boucle des scrutations entre debut et fin choisis' ; while (klu le kfin) do begin while (klu le kfin + 6) do begin ; l'ajout de 6 est destine a taiter le cas ou on veut isoler une ; seule scrutation (dans un fichier 10sec ou 128 sec), et entrer ; dans la boucle si une une frequence autre que la frequence 0 ; est selectionnee. ; Rem : pour une simulation la boucle est decrite une seule fois ; car on a pose plus haut kfin=1 et on a saute (jusqu'en a0:) ce ; qui pouvait le modifier. => n_image=1 pour une simulation. t0 = systime(1) ; temps en sec depuis le 1er jan 70. i_integ = 1 ; compteur d'integration de scrutations succes- ; sives (sera incremente jusqu'a fac_integ). b1: J_M = 0 ; fin du saut arriere pour integrer fac_integ ; scrutations pt = fltarr(2, 2, nbre_harm) ; initialisation de l'accumulateur RH_READCH, klu, sel_freq ; print, 'RH_DPATCHFITS_NRH, klu, kfin :', klu, kfin if klu eq -999L then goto, ecrit ; vers ecriture. if (i_simulation eq 0) then hmsc = bfi.h ch_format = "(i2.2)" ch_hms = string(format=ch_format, hmsc(0)) + ':' + $ string(format=ch_format, hmsc(1)) + ':' + $ string(format=ch_format, hmsc(2)) + ':' + $ string(format=ch_format, hmsc(3)) ; print, 'RH_DPATCHFITS_NRH, heure lue : ', ch_hms, $ ; ' num freq :',bfi.nof ; print, 'n_image =', byte(n_image) nf = bfi.nof ; numero de la frequence lue. ; print, klu, hmsc ; heure, min sec, 1/100 freq = float (entfi.frq(nf)) / 10 ch_format = "('RH_DPATCHFITS_NRH, heure lue : ', " + $ "i2, ':', i2, ':', i2, ' freq =', f6.1, ' MHz')" ; print, ' ' ; print, format=ch_format, hmsc(0), hmsc(1), hmsc(2), freq ; impressions souvent utile nsa = long (bfi.c) ; nbre d'images comprimees. ihms = transpose(temps)#hmsc ; heure en msec. ihms_lu = ihms ; pour controle local ; niveau de la boucle de lecture des scrutations successives. if (numrec(nf) eq 0L) then begin imstd(nf) = ihms itd(*,nf) = hmsc endif pt_lu = lufi.pt if (icorrec_phase eq 1) then RH_COREC_PHASE, nf, pt_lu, hmsc ; correction de l'erreur sur le phasage jusqu'au 4 nov 98. pt_lu = reform (pt_lu, 2 , 2 , entfi.nval / 4) ; pt_lu fltarr(4, 576 ou 648) puis fltarr(2, 2, 576 ou 648). ; Table des harmoniques : ; 1er indice: cos ou sin, ; 2eme ind: non polar ou polar, ; 3eme indice: balaie les 576 ou 648 harm. pt = pt + pt_lu / fac_integ if (i_integ eq 1) then ihms_1 = ihms if (i_integ eq fac_integ) then ihms_2 = ihms if (i_integ lt fac_integ) then begin i_integ = i_integ + 1 klu = klu + 1L ; numero de l'engrt suivant a lire, goto, b1 endif if (i_integ eq fac_integ) then begin ; print, 'fin d''integration (fac_integ =', fac_integ, ')' ihms = (ihms_1 + ihms_2) / 2 ; print, 'ihms_lu =', ihms_lu ; print, 'ihms =', ihms endif ; Debut du traitement de calcul d'image ; Conversion en canaux interferometriques de k_max et l_max (definis ; dans RH_DPNEW sous les noms de x_centrage et y_centrage) ; Cette conversion se fait ici puisqu'on dispose de l'heure hmsc. ; Elle n'est utile que pour un recentrage choisi (i_recentrage=2). if (i_recentrage eq 2) then begin x_centrage = x_centrage_0 y_centrage = y_centrage_0 if (version_solarsoft eq 0) then $ print, 'recentrage en coord helio (unites de 0.01 Rs) :', $ x_centrage, y_centrage ; Calcul des coordonnees du point de recentrage dans le repere ; (cosd *dH, ddelta), et transformation en radians ; rappel : p = (At, As) compte dans le sens trigo usuel. ; rappel : delta est defini et calcule plus haut. common MALAX, $ ; seul usage dans RH_MALC_IM_2D dans DPATCHFITS ij, im, ian, icorpoi, icorion, isour, $ etmer, edec, dew, new, hew, u, dns, hns, a1, nns, sinl, cosl, $ sinp, cosp, gdel, ghmer, ahmer, c, rsol, ahmersol, sind, cosd x_cen = x_centrage * cosp - y_centrage * sinp y_cen = x_centrage * sinp + y_centrage * cosp r_sol = !pi/180 / 60 * 16 x_cen = 0.01 * x_cen * r_sol y_cen = 0.01 * y_cen * r_sol ; Calcul du dephasage pour les bases dew et dns ; phi = 2*pi/lambda * base . d(vecteur unitaire de visee) base_ew = dew * [cos(u) * cos(hew), - cos(u) * sin(hew), sin(u)] base_ns = dns * [ -sinl * sin(hns), - sinl * cos(hns), cosl] ah = hmsc(0) + float(hmsc(1)) / 60 + float(hmsc(2)) / 3600 ang_hor = !pi/12 * (ah - ahmer) sinh = sin(ang_hor) cosh = cos(ang_hor) d_vec = [ cosh * x_cen - sinh * sind * y_cen, $ -sinh * x_cen - cosh * sind * y_cen, $ cosd * y_cen] lambda = c / freq d_phi_ew = 2*!pi / lambda * total(base_ew * d_vec) d_phi_ns = 2*!pi / lambda * total(base_ns * d_vec) x_cen_can = new / (2*!pi) * d_phi_ew ; canaux "interfero". y_cen_can = nns / (2*!pi) * d_phi_ns ; ----------------- x_centrage = x_cen_can y_centrage = y_cen_can endif if (version_solarsoft eq 0) then begin if (i_pave_plein eq 1) then begin print, ' Remplissage du pave central :' if (i_recentrage eq 0) then $ print, ' pas de recentrage' if (i_recentrage eq 1) then $ print, ' recentrage automatique sur maximum' if (i_recentrage eq 2) then $ print, ' recentrage en canaux interfero (de -npi/2 ' + $ 'a + npi/2 - 1)', byte(x_centrage), byte(y_centrage) endif endif ; Calcul de la maille elementaire interferometrique et de son angle solide. RH_MAILLE_3, npi, freq, hmsc, $ ; entree maille, domega ; sortie ; entree: npi nbre de pts sur l'image interferometrique. ; frequence (MHz), ; heure (h, m, s, c). ; common MALAX, contenant entre autres new (defini par ie0, ; ie1, ie2) et nns (toujours 64). ; sortie: maille = [x0, y0, x1, y1, x2, y2] contient les coord ; heliographiques en RS des points M0 (r=0, s=0), ; M1 (r=1, s=0), M2 (r=0, s=1) definissant la maille ; interferometrique MI. ; Ces coordonnees heliographiques tiennent compte des : ; . derives du soleil en angle horaire et declinaison, ; . rotation d'angle -p ou non selon que rot_p (de- ; fini en entree de rh_dpnew) vaut 1 ou 0. Ceci se ; fait par le biais de l'utilisation du common ; MALAX, qui est rempli par INIT_MALAX_2D, qui a ; rot_p dans sa liste d'entree. ; ; Un pavage de npi*npi MI recouvre le champ de l'image ; interferometrique calculee par fft en EW et NS. Ce ; champ peut etre double du champ interferometrique ; stricto sensu si ie0_ew = ie0_ns = 0. ; Rappel : npi est fige egal a 128 dans DPATCHFITS ; depuis jul 01. ; domega, angle solide sous-tendu par la maille (sterad). ; niveau de la boucle de lecture des scrutations successives. ; Calcul de la grille des coordonnees interferometriques correspondant a ; une grille heliographique carree et calcul du drapeau de limita- ; tion du champ "champ" dans le calcul de l'image. RH_GRILLE, npi, np, maille, n_period, champ_helio, freq, $; entree rg, sg, champ ; sortie. ; champ_helio (Rs) est la largeur totale du champ couvert par la ; grille heliographique a mailles carrees dans laquelle le soleil ; est rond. ; rg et sg sont fltarr(np, np) contenant les coord interferome- ; triques non entieres correspondant aux noeuds de la grille ; des coordonnees heliographiques `a maille carree sur laquelle ; on calcule l'image du soleil. ; On tient toujours compte des derives en angle horaire et decli- ; naison du soleil, et de l'angle p seulement si rot_p=1 en ; entree de RH_DPNEW. ; "champ" en sortie : drapeau (n*n) qui vaut 1 pour les couples ; (rg, sg) interieurs a [0, n], c'est a dire au champ interfe- ; rometrique tenant compte des harmoniques de E0. ; Rappel: ; Les indices interferometriques r et s varient sur [0, npi] sur ; un champ interferometrique. Si une image periodisee ils ; varient sur [0, n_period*npi]. ; rg et sg peuvent varier sur un domaine plus etroit ou plus lar- ; ge selon que le champ champ_helio choisi par l'utilisateur ; est plus petit ou plus grand que le champ interferometrique ; eventuellement periodise. ; champ est un intarr(np, np), comme rg et sg, et qui vaut zero la ; ou rg et sg correspondent a des points exterieurs au champ ; interferometrique eventuellement periodise. ; Avant le 8 sept 99 "champ" valait aussi zero pour r < rmax du ; centre solaire (rmax defini dans RH_GRILLE). C'est une de- ; louiserie abandonnee. ; Visualisation des tableaux rg et sg. Voir dans RH_GRILLE ; niveau de la boucle de lecture des scrutations successives. ; Calcul du tableau dom_c_i des indices des points de l'image interferometrique ; (domaine qui sera effectivement pris en compte dans le clean et ; etendu a toute l'image si mode_clean=1) ; dom_c_i comme "domaine C interferometrique". RH_CALC_DOM_C, $ mode_clean, x_dom, y_dom, dr_dom, $ ; entree npi, maille, $ ; entree dom_c_i ; sortie ; Definition de dom_c_i suivant les valeure de mode_clean : ; . 1 (clean usuel) : totalite de l'image interferometrique ; (npi*npi). ; . 2 domaine de non restitution des composantes, recherchees ; sur toute l'image. ; . 3 domaine de suppression des composantes (debarasser l'ima- ; d'un centre intense). ; . 4 domaine de nettoyage (amelioration de l'image d'un centre ; intense. icon = 0 if (icon eq 1) then begin toto = intarr(npi, npi) toto(dom_c_i) = 1 ; Trace des axes et du cadre toto_512 = congrid(toto, 512, 512, cubic=-0.5) toto_512(*, 510) = 1 toto_512(510, *) = 1 toto_512(*, 256) = 1 toto_512(256, *) = 1 window, 0 & wset, 0 shade_surf, toto window, 1 & wset, 1 tvscl, toto_512 stop endif ; niveau de la boucle de lecture des scrutations successives. ; Calcul du tableau dom_a_i des indices des points de l'image interferometrique ; destines a visualiser le contour de dom_c defini dans DPNEW ; (ce contour est toujours figure sur les traces de controle de ; RH_CALC_IMAGE_HELIO par commodite, meme si mode_clean=1 et si ; dom_c est ensuite etendu a toute l'image) ; dom_a_i comme "domaine annulaire interferometrique". mode = 2 RH_CALC_DOM_C, $ mode, x_dom, y_dom, dr_dom + 0.02, $ ; entree npi, maille, $ ; entree dom_ext ; sortie RH_CALC_DOM_C, $ mode, x_dom, y_dom, dr_dom - 0.02, $ ; entree npi, maille, $ ; entree dom_int ; sortie flag_e = intarr(npi, npi) flag_i = intarr(npi, npi) flag_e(dom_ext) = 1 flag_i(dom_int) = 1 dom_a_i = where(flag_e ne flag_i) icon = 0 if (icon eq 1) then begin toto = intarr(npi, npi) toto(dom_a_i) = 1 window, 0 tvscl, congrid(toto, 512, 512, cubic=-0.5) stop endif ; Calcul du tableau dom_a_h des indices des points de l'espace annulaire circu- ; laire) entourant dom_c sur l'image heliographique (contour tou- ; jours figur'e par commodite sur les traces, meme pour ; mode_clean=1) ; dom_a_h comme "domaine annulaire heliograpique". d_o = shift(dist(np), -np/2, -np/2) * champ_helio / np ; d_o dist a l'origine (Rs). Rappel : champ_helio est en Rs. conv = float(np) / float(champ_helio) ; conversion RS -> pixels. x_pix = fix(conv * x_dom + 0.5) ; absc centre CME en pixels. y_pix = fix(conv * y_dom + 0.5) ; ord -------------------- d_c = shift(d_o, x_pix, y_pix) ; distance au centre de dom_c ; en RS mais x_pix et y_pix ; sont en pixels. dom_a_h = where((d_c le dr_dom + 0.03) and (d_c ge dr_dom - 0.03)) ; Ici figuraient quelques delouiseries sur la fumeuse variable "champ" ; ainsi que des commentaires tentant de rendre les choses comprehensi- ; bles (entreprise a priori vouee a l'echec). Le tout est transfere a ; la cave (plus bas). ; Rem : on est toujours dans la boucle des scrutations successives. ; niveau de la boucle de lecture des scrutations successives. ; Aiguillage eventuel vers polar seul if (i_polar eq 2) then goto, a1 ; Traitement non polar ; Remplissage du tableau des harmoniques non polar rang'es format Nancay. harm_N_np = complex(pt(0, 0, *), pt(1, 0, *)) ; le 2eme indice indique non polar ou polar. harm_N_np = reform(harm_N_np) ; pour supprimer 2 indices degeneres. ; pt est un reel, tableau des harmoniques. 1er indice: cos ou sin, ; 2eme ind: non polar ou polar, 3eme indice: balaie les 576 ou 648 ; harm. harm_N_np est donc un complexarr (576 ou 648). ; Correction des irregularites d'implantation d'antennes if (i_corr_pos_ant eq 1) then $ CORR_POS_ANT, $ ; Dans ~/calib/rh_corr_antennes.pro n_image, $ ; entree i_source_cal, h_cal, $ ; ------ date_obs, dec, hmer, hmsc, $ ; ------ correl, freq, $ ; ------ harm_N_np ; entree et sortie ; . i_source_cal identifie la source de calibration, dont le ; l'heure meridien et la declinaison sont definies dans ; CORR_POS_ANT (precision non necessaire). ; . date_obs permet de savoir s'il faut ou non retrancher les ; "corrections" Delouis et si observation temps reel et ; calibration ont utilise des positions corrigees (et exac- ; tes) des antennes ou des positions approchees deduites ; des reseaux moyens, ce qui sera le cas au moins jusque ; vers nov 2003. Claude previendra du changement. ; . dec et hmer pour l'objet observ'e, en TU pour le Soleil et ; en TS pour les sources de controle. ; . hmsc : 2 premiers elements seuls utilises dans CORR_POS_ANT ; Rappel : on a decide, par simplicite, que le Temps Reel et ; la calibration (ou la recalibration) utilisaient les ; memes valeurs fausses des positions d'antennes. Ceci ; a l'avantage que la correction de position d'antennes ; est localisee uniquement dans CORR_POS_ANT et que la ; formule du dephasage de correction est unique. ; On espere que des valeurs fiables et definitives des ; positions d'antennes seront disponibles debut 2004. ; A ce moment onn les introduira dans le Temps Reel et la ; calibration et CORR_POS_ANT se cout-circuitera de lui ; meme au vu de la date d'observation. ; Rem : COR_POS_ANT corrige une deformation du lobe variable ; avec le temps. Le lobe resultant est donc presque ; invariable mais il peut etre tres moche. ; Correction de la rotation des dipoles et/ou du depointage des antennes ; lunettes (utile pour les sources de controle) ; date_obs suffit pour savoir comment corriger. Le cas ; de l'examen de sources de calibration n'est pas prevu au 11 avr ; 03. ; if (i_con_cal eq 1) then begin ; seul cas des sources de controle. type_obs = 1 & type_calib = 1 ; bidons. CORR_ROT_DEP_NP, $ ; fichier RH_CORR_ANTENNES n_image, $ ; entree. date_obs, type_obs, type_calib, $ ; ------ rotew, rotns, ah_dep, delta, freq, $ ; ------ (rads, MHz). harm_N_np, $ ; ------ et sortie. a_H2, a_NS0 ; sortie (radians), endif ; Controle des harmoniques NS de NS0_ew et NS0_ns icon = 0 if (icon eq 1) then begin ch_polar = 'non polar : ' DP_CON_1, freq, ch_polar, harm_N_np endif ; niveau de la boucle de lecture des scrutations successives. ; Remplacement eventuel des donnees reelles par des donnees simulees ; Rappel : la simulation concerne ici uniquement le cas solaire ; puisque la calibration est entierement faite dans RH_DP_CALIB. if (i_simulation eq 0) then begin objet_simule = 0 ; au lieu de fltarr(npi, npi) pour economie. tf_objet_simule = 0 ; ----------------------------------------- ; objet_simule et tf_objet_simule sont bidons pour figurer ; dans la liste de RH_CALC_IMAGE_HELIO mais ne seront pas ; utilises. endif if (i_simulation gt 0) then begin i_choix_sim = i_simulation ; vestige d'avant le 21 fev 06, quand ; le choix de l'objet simul'e n'etait ; pas remont'e dans rh_dpnew. Rappel : ; . 1 pour genre Cygne, et la simula- ; tion est dans l'espace reel. ; . 2 pour soleil +- complexe, et la ; simulation est dans l'espace ; interferometrique. ; print, 'RH_DPATCHFITS_NRH : appel SIMUL_HARM_N donnees Soleil' ; Rappel : apres le 18 mars 03 le calcul de la calibration (ou ; de la recalibration) elle meme est elimin'e. Il ne reste ; que l'usage eventuel des fichiers coeff_corr (dont le nom ; est ici contenu dans la variable ch_cor) produits par ; RH_DP_CALIB . i_calib = 1 ; indique (pour une simulation genre Cygne) que ; l'objet est simul'e dans l'espace reel (et ; non "interferometrique") ; pour une simulation genre soleil (i_choix_sim ; = 2) i_calib est inutilis'e dans ; SIMUL_HARM_N SIMUL_HARM_N, $ ; dans RH_SIMULATION. ; entree i_choix_sim, i_calib, $ freq, t_moy, correl, $ ; i_pave_plein, i_sys_lin, n_bord, $ ; controles seuls ; sortie gain_ant, phase_ant, $ ; objet_simule, pix_obj, $ ; objet, pixel objet. tf_objet_simule, mat_u, mat_v, $ ; visib sur la grille ; RH, coord des pts ; de la grille. harm_N_np ; i_sys_lin est pass'e a RH_MALC_IM_2D pour controles. ; phase_ant est en radians ; mat_u et _v coord des pts de la grille (matrices 128*128) ; tf_objet_simule est centree en milieu de tableau. ; rappels : ; . la simulation des donnees pour une calibration est ; faite apres integration, apres endwhile. ; . n_bord est necessaire a MALC_IM_2D, pour controles. ; . i_calib (ici 0 ou 1) est indifferent car on appelera ; SIMUL_VISIB_2. Il n'y a pas de confusion avec le ; i_calib de RH_DP_CALIB (qui appelle SIMUL_VISIB_1). endif ; fin du if de simulation non polar. ; niveau de la boucle de lecture des scrutations successives. ; Controle de l'objet et de sa TF icon = 0 if (icon eq 1) then begin help, i_simulation, objet_simule, pix_obj, tf_objet_simule window, 0 loadct, 0 objet_512 = congrid(objet_simule, 512, 512, cubic=-0.5) tvscl, objet_512 window, 1 shade_surf, abs(tf_objet_simule) stop endif ; Controle des harm de NS8 (2 fev 2006) icon = 0 if (icon eq 1) then begin harm_ns8 = harm_N_np (540 : 557) window, 0 plot, abs(harm_ns8) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'abs(harm avec NS8)' ch_format = "(18i4)" print, '' print, 'correl (0, 540:557) :' print, format=ch_format, correl(0, 540:557) print, '' print, 'correl (1, 540:557) :' print, format=ch_format, correl(1, 540:557) stop ; conclusion : OK c'est bien symetrique par rapport a NS8. endif ; Correction eventuelle de l'image non polar a l'aide d'une recalibration ; a posteriori (directe ou auto) faite avec un Cygne a une date ; proche de celle de l'observation, et avec les memes tables de ; correction ; Rappel : la recalibration (comme la calibration) est faite avec ; les memes positions d'antennes (fausses ou justes) que ; le Temps Reel. Cela simplifie la gestion des correc- ; tions, faites en un seul endroit (CORR_POS_ANT). if (i_recalib eq 1) then begin i_np_p = 0 ; 0 pour polar, 1 polar. RH_RECALIBRATION, $ ; fichier ~/calib/rh_recalibration.pro n_image, $ ; entree i_np_p, freq, fich_cor, $ ; ------ correl, $ ; ------ harm_N_np ; ------ et sortie endif ; Controle des harm de NS8 (2 fev 2006) ; a permis de troiuver une faute dans la correction des harm de E0 ; dans CORR_GAIN_ANT. icon = 0 if (icon eq 1) then begin harm_ns8 = harm_N_np (540 : 557) window, 0 plot, abs(harm_ns8) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'abs(harm avec NS8)' ch_format = "(18i4)" print, '' print, 'correl (0, 540:557) :' print, format=ch_format, correl(0, 540:557) print, '' print, 'correl (1, 540:557) :' print, format=ch_format, correl(1, 540:557) stop ; conclusion : OK c'est bien symetrique par rapport a NS8. endif ; Correction eventuelle de certaines antennes (gains entres manuellement) ; Il n'y a pas de "if" car on entre le drapeau g_flag. phi_ew_ns_np = phi_ew_ns(0) & phi_w_e_np = phi_w_e(0) g_amp_np = g_amp(0) & g_phi_np = g_phi(0) CORR_GAIN_ANT, $ ; fichier RH_CORR_ANTENNES.PRO n_image, $ ; entree. phi_ew_ns_np, phi_w_e_np, $ ; ------ g_flag, g_amp_np, g_phi_np, $ ; ------ correl, $ ; ------ harm_N_np ; entree et sortie ; niveau non polar de la boucle de lecture des scrutations successives. ; Controle des harmoniques NS avec NSO_ew et NS0_ns (22 fev 2006) ; NS0 ---- NS10 a NS23 et NS45 558 a 572 ; --- ---- NS7 a NS9 573 a 575 icon = 0 if (icon eq 1) then begin h_NS_NS0_ew = complexarr(25) h_NS_NS0_ns = complexarr(25) for ii = 1, 23 do h_NS_NS0_ew(ii) = harm_N_np (9 + 18*ii) h_NS_NS0_ns (10) = harm_N_np (558 : 572) h_NS_NS0_ns ( 7) = harm_N_np (573 : 575) amax = max(abs(h_NS_NS0_ew)) > max(abs(h_NS_NS0_ns)) loadct, 0 window, 1 plot, abs(h_NS_NS0_ew), xrange=[-1, 25], xstyle=1, $ yrange=[0, 1.2 * amax], ystyle=1 oplot, 0.98 * abs(h_NS_NS0_ns), linestyle=1 xyouts, 0.5, 0/97, /normal, alignment=0.5, charsize=1.3, $ 'abs des harm NS avec NS0_ew ( ___ ) et ns0_ns (...) stop endif ; Calcul de l'image non polar et normalisation ; Le temps cpu utilis'e depuis le debut de la boucle sur les scru- ; tations est environ 0.162 sec sur msopz ; print, 'RH_DPATCHFITS_NRH : appel RH_CALC_IMAGE_HELIO' ; t1 = systime(1) ; temps en sec depuis le 1er jan 70. ; ch_1 = "(' Temps depuis debut boucle sur les scrutations', " ; ch_2 = "' sec')" ; ch_format = ch_1 + " f10.3," + ch_2 ; print, format=ch_format, t1 - t0 ch_pol = 'non polar' anorm_1_s_np = 1 ; bidon car sera ecras'e. anorm_1_p_np = 1 ; bidon car sera ecras'e. RH_CALC_IMAGE_HELIO, $ ; Entree ; Caracterisation image jul_dat, correl, freq, n_image, npi, n_bord, $ i_polar, ch_pol, anorm_1_s_np, anorm_1_p_np, $ n_period, mode_interpol, np, champ_helio, t_b, $ i_affich_lobe, $ ; 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, $ ; Complement du pave central (harm multiples impairs de 50m) i_pave_plein, i_recentrage, x_centrage, y_centrage, $ ; Donnees pour une image i_simulation, tf_objet_simule, objet_simule, $ harm_N_np, rg, sg, domega, champ, $ ; Clean mode_clean, dom_c_i, dom_a_i, dom_a_h, $ i_ech_min, i_ech_max , a_in, a_ex, $ a_red_fil, itermax, gain_iter, jmax, crit_arret, $ i_seuil, coeff_seuil, frac_flux, fac_crois_res, $ i_test_arret, i_stop_image, i_ecrete, frac_ecr, $ tete_fich, $ i_lissage, poids_c, larg_2, $ ; Divers i_sys_lin, $ ; Sortie: x_centrage_np, y_centrage_np, $ l_th_np, im_2d_np, f_total_np, f_compact_np ; En sortie : ; im_2d_np image ronde non polar, corrigee : ; . des derives du soleil en angle horaire et declinaison, ; . de la rotation d'angle -p seulement si rot_p=1 en ; entree de RH_DPNEW. ; normalisation telle que : ; . somme des points = flux en sfu si t_b=1, ; . brillance en kelvins si t_b=2. ; f_total est le flux total en sfu. C'est la somme des points de ; l'image interferometrique sale obtenue par fft dans ; RH_MALC_IM_2D (la base etant ajustee de facon que la ; moyenne de l'image soit nulle sur les bandes N et S ; de l'image. On verifie que la somme des points de ; l'image interferometrique propre est tres voisine. ; La transformation en image heliographique est faite en- ; suite : ; . si t_b=1, l'image helio est normalisee de facon que ; la somme de ses points donne aussi flux_total. ; . si t_b=2, l'image helio est normalisee de facon ; qu'elle soit en Kelvins. ; Au 28 nov 2005 f_total_np et f_total_p sont inuti- ; lises dans la suite. ; f_compact est le flux de la source compactes, obtenu par extapola- ; tion quadratique de la visibilite dans RH_MALC_IM_2D. ; x_centrage_np est l'abscisse de recentrage automatique non polar ; (canaux inteferometriques) pour passage au cas ; polar, pour lequel le recentrage est considere'e ; comme moins fiable. ; Rappels : ; . rg et sg sont les coord interfero des noeud d'une grille ; helio a mailles carrees. Dans leur calcul on a tenu compte: ; . des derives du soleil em angle horaire et declinaison, ; . de la rotation d'angle -p seulement si rot_p=1 en ; entree de RH_DPNEW. ; . tf_objet_simule est centree en milieu de tableau. ; . le calcul du tableau dom_c_i est fait juste apres GRILLE. ; niveau non polar de la boucle de lecture des scrutations successives. icon = 0 if (icon eq 1) then begin max_im_helio = max(im_2d_np) ch_max = string(format="(g8.2)", max_im_helio) if (t_b eq 2) then ch_max = ch_max + ' K' chsiz = 1.5 i_sh_tv = 2 ; 1 shade_surf, 2 tvscl im_visu = congrid(im_2d_np, 512, 512, cubic=-0.5) l_Th_visu = congrid(l_th_np , 512, 512, cubic=-0.5) if (i_sh_tv eq 1) then loadct, 0 if (i_sh_tv eq 2) then loadct, 13 ; 13 eainbow window, 0 & wset, 0 if (i_sh_tv eq 1) then shade_surf, l_th_visu if (i_sh_tv eq 2) then tvscl , l_th_visu xyouts, /normal, 0.5, 0.95, alignment=0.5, charsize=chsiz, $ 'lobe theorique "interferometrique" sortie calc_image_helio' window, 1 & wset, 1 if (i_sh_tv eq 1) then shade_surf, im_visu if (i_sh_tv eq 2) then tvscl , im_visu xyouts, /normal, 0.5, 0.95, alignment=0.5, charsize=chsiz, $ 'image heliographique sortie calc_image_helio, max : ' + $ ch_max stop endif max_np = max(im_2d_np, ind_np, min=min_np) ; rappel: n_image est le numero de l'image traitee. ; Rem: le 11 jan 2 000 il est verifie que l'usage de congrid jusqu'a ; la sortie de CALC_IMAGE_2D est limite a des verification. ; Rappel: quand on passe de 128 a 512 pts congrid calcule l'image ; sur un champ un peu plus grand, comme INTER_FFT. Les pts ; de cette partie supplementaire du champ sont toutefois ; faux par rapport a ceux de RH_INTER_FFT, mais pour le reste ; c'est tres voisin. ; print, 'RH_DPATCHFITS_NRH : heure scrut precedemment traitee : ', $ ; ch_hms ; Fin de calcul de l'image non polar et normalisation ; Fin de traitement non polar ; niveau non polar de la boucle de lecture des scrutations successives. a1: J_M = 0 ; fin du saut de l'image non polar si i_polar=2 . ; Debut de traitement de l'image polar if ((i_polar eq 1) or (i_polar eq 2)) then begin harm_N_p = complex(pt(0, 1, *), pt(1, 1, *)) ; le 2eme indice indique non polar ou polar. harm_N_p = reform(harm_N_p) ; pour supprimer 2 indices degeneres. if keyword_set(selfga_p) then harm_N_p = harm_N_p * selfga_p ; niveau polar dans boucle de lecture des scrutations successives. ; Correction des irregularites d'implantation d'antennes if (i_corr_pos_ant eq 1) then $ CORR_POS_ANT, $ ; Dans fichier RH_CORR_ANTENNES.PRO n_image, $ ; entree i_source_cal, h_cal, $ ; entree date_obs, dec, hmer, hmsc, $ ; ------ correl, freq, $ ; ------ harm_N_p ; entree et sortie ; Commentaires sur les arguments : voir appel non polar ; Correction des harmoniques dans le cas de l'examen des sources de ; controle ; Dans ce cas date_obs suffit pour savoir comment corriger. Le ; cas de l'examen de sources de calibration n'est pas prevu ; au 11 avr 03. if (i_con_cal eq 1) then begin type_obs = 1 & type_calib = 1 CORR_ROT_DEP_P, $ n_image, $ ; entree. date_obs, type_obs, type_calib, $ ; ------ rotew, rotns, ah_dep, delta, freq, $ ; ------ (rads, MHz). harm_N_p, $ ; ------ et sortie. a_H2, a_NS0 ; sortie (radians), endif ; niveau polar dans boucle de lecture des scrutations successives. ; Controle des harmoniques NS de NS0_ew et NS0_ns icon = 0 if (icon eq 1) then begin ch_polar = 'polar : ' DP_CON_1, freq, ch_polar, harm_N_p endif ; niveau de la boucle de lecture des scrutations successives. ; Remplacement eventuel des donnees polar par des donnees simulees ; La simulation concerne ici uniquement les donnees solaires. ; Dans le cas de la calibration elle est faite plus bas. if (i_simulation eq 0) then begin objet_simule = 0 ; au lieu de fltarr(npi, npi). Economie. tf_objet_simule = 0 ; ------------------------------------- endif if (i_simulation gt 0) then begin i_choix_sim = i_simulation ; vestige d'avant le 21 fev 06, quand le choix de l'objet ; simul'e n'etait pas remont'e dans rh_dpnew. Rappel : ; . 1 pour genre Cygne, . 2 pour soleil +- complexe. i_calib = 1 ; bidon comme pour non polar avec ; i_choix_sim = 2 SIMUL_HARM_N, $ ; simul polar soleil complexe. ; entree i_choix_sim, i_calib, freq, t_moy, correl, $ i_pave_plein, i_sys_lin, n_bord, $; controles ; sortie gain_ant, phase_ant, $ objet_simule, tf_objet_simule, harm_N_p ; i_sys_lin est pass'e a RH_MALC_IM_2D pour controles. ; phase_ant est en radians ; tf_objet_simule est centree en milieu de tableau. ; rappels : ; . la simulation des donnees pour une calibration est ; faite apres integration, apres endwhile. ; . n_bord est necessaire a RH_MALC_IM_2D, pour controles. ; . i_calib (ici 0 ou 1) est indifferent car on appelera ; SIMUL_VISIB_2. Pas de confusion avec le i_calib de ; RH_DP_CALIB (qui appelle SIMUL_VISIB_1). endif ; fin du if de simulation polar. ; niveau polar dans boucle de lecture des scrutations successives. ; Correction eventuelle a l'aide d'un calibration de type auto if (i_recalib eq 1) then begin i_np_p = 1 ; 0 pour polar, 1 polar. RH_RECALIBRATION, $ ; fichier RH_RECALIBRATION.PRO n_image, $ ; entree. i_np_p, freq, ch_cor, $ ; ------ correl, $ ; ------ harm_N_p ; ------ et sortie endif ; niveau polar dans boucle de lecture des scrutations successives. ; Correction eventuelle de certaines antennes ; Il n'y a pas de "if" car on entre le drapeau g_flag. phi_ew_ns_p = phi_ew_ns(0) & phi_w_e_p = phi_w_e(0) g_amp_p = g_amp(0) & g_phi_p = g_phi(0) CORR_GAIN_ANT, $ ; fichier RH_CORR_ANTENNES.PRO n_image, $ ; entree. phi_ew_ns_p, phi_w_e_p, $ ; entree. g_flag, g_amp_p, g_phi_p, $ ; ------ correl, $ ; ------ harm_N_p ; entree et sortie. ; Calcul de l'image polar et normalisation ch_pol = 'polar' if (i_polar eq 1) then begin ; polar apres non polar if (i_pave_plein eq 0) then begin i_recentr_pol = 0 ; bidon x_centrage = 0 ; ---- y_centrage = 0 ; ---- endif if (i_pave_plein eq 1) then begin i_recentr_pol = 2 ; meme centrage que np x_centrage = x_centrage_np y_centrage = y_centrage_np endif endif if (i_polar eq 2) then begin ; cas polar seule. i_recentr_pol = i_recentrage ; on ne fait rien : le recentrage automatique le reste, et ; s'il y a des problemes de recentrage sur l'alias plutot ; que sur la region active, on recommence en prenant l'op- ; tion du recentrage choisi. endif RH_CALC_IMAGE_HELIO, $ ; Entree ; Caracterisation image jul_dat, correl, freq, n_image, npi, n_bord, $ i_polar, ch_pol, anorm_1_s_np, anorm_1_p_np, $ n_period, mode_interpol, np, champ_helio, t_b, $ i_affich_lobe, $ ; Suppression (eventuelle) harmoniques de certaines antennes ie0_ew, ie0_ns, ie1_ew, ie2_ew, ie2_ns, $ ins45_ew, ins45_ns, i_redond, gap_ns0, sel_ant, $ ; Complement du pave central (harm multiples impairs de 50m) i_pave_plein, i_recentr_pol, x_centrage, y_centrage, $ ; Donnees pour une image i_simulation, tf_objet_simule, objet_simule, $ harm_N_p, rg, sg, domega, champ, $ ; Clean mode_clean, dom_c_i, dom_a_i, dom_a_h, $ i_ech_min, i_ech_max , a_in, a_ex, $ a_red_fil, itermax, gain_iter, jmax, crit_arret,$ i_seuil, coeff_seuil, frac_flux, fac_crois_res, $ i_test_arret, i_stop_image, i_ecrete, frac_ecr, $ tete_fich, $ i_lissage, poids_c, larg_2, $ ; Divers i_sys_lin, $ ; Sortie: x_centrage_p, y_centrage_p, $ l_th_p, im_2d_p, f_total_p, f_compact_p ; im_2d_p image ronde et normalisee. ; Rappel : tf_objet_simule est centree en milieu de tableau. max_p = max(im_2d_p, ind_p, min=min_p) ; print, 'RH_DPATCHFITS_NRH (apres RH_CALC_IMAGE_HELIO' + $ ; ' polar) heure scrut traitee : ', ch_hms endif ; Fin de traitement de l'image polar. ; niveau dans boucle de lecture dscrutations successives (np et/ou p). ; Rafraichissement des maxima et minima (non polar ou polar) mmax et mmin ; pour la journee. Rappel : nf numero de frequence. if ((i_polar eq 0) or (i_polar eq 1)) then begin ; les max et min sont les max et min non polar. mmax(nf) = max([max_np, mmax(nf)]) mmin(nf) = min([min_np, mmin(nf)]) endif if (i_polar eq 2) then begin ; seul le polar est enregistr'e mmax(nf) = max ([max_np, max_p, mmax(nf)]) mmin(nf) = min ([min_np, min_p, mmin(nf)]) endif ; Comptage des enregitrements (un seul enreg fits contient les images ; non polar et polar pour une seule frequence et pour une heure ; donnee). C'est la valeur du numero de frequence nf qui aiguille ; l'ecriture vers le fichier correspondant a une frequence donnee. numrec(nf) = numrec(nf) + 1L ; numero de l'image a l'ecriture.; if (version_solarsoft eq 1) then begin ; version SolarSoft. retchar = string([13B]) ch_formatssw = "($,' Ecrit. num. ', i6, 4x,' Frq.',i2,2x, " + $ "i2.2, ':',i2.2, ':', i2.2, ':', i2.2, " + $ "' klu =', i6, ' kfin =', i6,a1)" if (numrec(nf) mod 10) eq 0 then print, format=ch_formatssw, $ numrec(nf), bfi.nof, bfi.h, klu, kfin, retchar end else begin ; version de Claude print, '' ch_format = "('Ecriture image(s) fits numero ', i6, ' :')" ; print, format=ch_format, numrec(nf) ch_format_1 = "(4x, 'derniere image lue = ', i2.2, " + $ "':', i2.2, ':', i2.2, ':', i2.2, ' klu =', i6, " + $ "' kfin =', i6)" ; print, format=ch_format_1, bfi.h, klu, kfin endelse if (i_simulation gt 0) then goto, a4 ; niveau dans boucle de lecture dscrutations successives (np et/ou p). ; Ecriture dans les fichiers mono-frequence ; on ecrit l'oppos'e im_2d_p2 de l'image polar im_2d_p si ; i_mul_pol = -1 (pour faciliter la relecture de fichiers fits ; polar quand la polar est negative). if (i_polar gt 0) then im_2d_p2 = i_mul_pol * im_2d_p ; Controle des images qu'on va ecrire icon = 0 if (icon eq 1) then begin window, 0, xsize=512, ysize=512 loadct, 13 tvscl, congrid (im_2d_np, 512, 512, cubic=-0.5) loadct, 0 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'image non polar' window, 1, xsize=512, ysize=512 loadct, 13 tvscl, congrid (im_2d_p2, 512, 512, cubic=-0.5) if (i_mul_pol eq 1) then ch_lab = 'non changee de signe' if (i_mul_pol eq -1) then ch_lab = 'changee de signe' loadct, 0 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'image polar ' + ch_lab max_im_np = max (im_2d_np) max_im_p = max (abs(im_2d_p)) print, 'controle RH_DPATCHFITS_NRH : ' print, ' maxima des valeurs absolues d''images helio (K) :' print, format="(4x, 'non polar :', e9.1, 6x, 'polar :', e9.1)", $ max_im_np, max_im_p stop ; icici endif fxbwrite, iunit(nf), ihms, 1, numrec(nf) fxbwrite, iunit(nf), numim(nf), 2, numrec(nf) if ((i_polar eq 0) or (i_polar eq 1)) then $ fxbwrite, iunit(nf), im_2d_np, 3, numrec(nf) if ((i_polar eq 1) or (i_polar eq 2)) then $ fxbwrite, iunit(nf), im_2d_p2, 4, numrec(nf) imstf(nf) = ihms itf(*, nf) = hmsc if (version_solarsoft eq 1) then goto, a6 ; version solarsoft ; Impression de recapitulation hmsc_fits = intarr(4) hmsc_fits (0) = ihms / 3600000L hmsc_fits (1) = (ihms - 3600000L * hmsc_fits (0)) / 60000L hmsc_fits (2) = (ihms - 3600000L * hmsc_fits (0) - $ 60000L * hmsc_fits (1)) / 1000L hmsc_fits (3) = (ihms - 3600000L * hmsc_fits (0) - $ 60000L * hmsc_fits (1) - $ 1000L * hmsc_fits (2)) / 10 ch_format = "('Ecriture image fits : heure = ', i2.2, ':', i2.2, " + $ "':', i2.2, ':', i2.2, 4x, '(numero ', i6, ')')" print, format=ch_format, hmsc_fits, numrec(nf) print, ' flux (sfu) total compact max (K)' + $ ' min (K)' if ((i_polar eq 0) or (i_polar eq 1)) then begin max_np = max((im_2d_np), min = min_np) ch_format = "(2x, 'non polar :', 10x, f7.1, 2x, f7.1, 4x, " + $ "e8.1, 4x, e8.1)" print, format=ch_format, f_total_np, f_compact_np, max_np, min_np endif if ((i_polar eq 1) or (i_polar eq 2)) then begin max_abs_p = max(abs(im_2d_p)) max_p = max(im_2d_p, min = min_p) if (max_p eq max_abs_p) then i_sign_pol = 1 else i_sign_pol = -1 if (i_sign_pol eq 1) then max_p = max (im_2d_p, min = min_p) if (i_sign_pol eq -1) then max_p = min (im_2d_p, max = min_p) ; si la polar est negative on imprime - maximum de la valeur ; absolue, c'est a dire le minimum algebrique. ch_format = "(2x, ' polar :', 10x, f7.1, 2x, f7.1, 4x, " + $ "e8.1, 4x, e8.1)" print, format=ch_format, f_total_p, f_compact_p, max_p, min_p endif print, '' print, '' a6: J_M = 0 ; fin du saut des impressions de controle. ; Calcul du numero suivant d'image a ecrire (ci-dessus, au tour suivant) nsa = long(bfi.c) numim(nf) = numim(nf) + nsa + 1L ; print, 'RH_DPATCHFITS_NRH : numero suivant d'image a ecrire =', numim a4: j_M = 0 ; fin du saut de calcul de rafraichissement du maximum ; et d'ecriture dans le cas de la calibration. klu = klu + 1L ; numero de l'engrt suivant a lire, n_image = n_image + 1 ; nbre d'image suivante, toutes frequences ; jointes. if ((i_simulation eq 0) and (n_image gt numtot)) then goto, ecrit ; vers ecriture (qui sera sautee pour des donnees reelles ou une ; calibration. endwhile ; Fin de boucle sur les acquisitions et les frequences if (i_simulation gt 0) then goto, a5 ; Fermeture des fichiers d'image separes en frequence, fermeture fichiers lus. ecrit: J_M = 0 ; Ecriture des fichiers d'image separes en frequence n_image = n_image - 1 ; nbre reel d'image (dernier +1 en trop). RH_CLOSE_HEADER, $ ; fermeture des fichiers d'ecriture. iunit, dirdsmf, entfi.nf, sel_freq, $ itd, itf, numrec, np, $ entfi.dat, imstd, imstf, mmin, mmax ; print,' sortie de rh_close_header',entfi.klumax a5: J_M = 0 RH_CLOSE ; ferme le fichier de lecture. print, 'sortie de RH_DPATCHFITS_NRH, fin de traitement et ecriture fits' if (version_solarsoft eq 0) then stop ; pour ne pas perdre les variables. end ; fin de RH_DPATCHFITS_NRH.PRO ;______________________________________________________________________________ ; GALERIE ARCHEOLOGIQUE CONSACREE AUX DELOUISERIES ; Delouiseries (interet purement historique) ; Calcul des tableaux des abscisses et ordonnees normalisees (comprises entre ; -1 et 1-1/64) des points de l'image. Ils ne servent qu'au calcul de ; limitation du champ facon J-M (cas i_lim_champ=2). ; xn = fltarr(np, np) ; grille des abscisses des points de l'image. ; for i = 0, np - 1 do xn(i, *) = (i - 0.5*np) / (0.5*np) ; xn varie de -1 a 1 - 1/64 ; yn = transpose (xn) ; Remarques: ; - Le point (0,0) de l'image est en bas a gauche (coin SE du soleil), ; alors qu'il est en haut a gauche pour print, tab. Ainsi x croit ; dans une ligne et y croit avec le numero de la ligne. ; - Les tableaux x(n,n) et y(n,n) des coordonnees (normalisees entre ; -1 et +1) des points de l'image (dans laquelle le soleil est rond) ; ont resp. leurs lignes et leurs colonnes identiques. Ils se dedui- ; sent l'un de l'autre par transposition. Ces points sont par choix ; sur une grille carree. Les valeurs de l'image en ces points sont ; obtenues plus bas par une interpolation de l'image ou le soleil ; n'est pas rond (elle meme obtenue par fft), aux points de coord ; interferometriques ("canaux") non entieres et calculees par GRILLE ; (rij et sij). ; Visualisation du drapeau de permission de calcul d'image. ; icon = 0 ; if (icon eq 1) then begin ; help, champ ; shade_surf, congrid(champ, 512, 512, cubic=-0.5) ; stop ; endif ; Fin de controle. ; Variable definie pour contouner une delouiserie (maintenant au placard) ; i_lim_champ = 1 ; 1 dist helio du centre du disque < rmax ; (rmax = 1.7 RS dans GRILLE). Le travail ; est deja fait dans GRILLE. ; 2 facon J-M (contestable). ; Commentaires judicieux ; if (i_lim_champ eq 2) then begin ; debut facon J-M. ; champ = where ( ((xn^2 + yn^2)^0.5) lt (1.7*164/freq) ) ; ; Rappel: xn et yn varient de -1 a 1-1/64. La limitation du ; ; champ est donc a 3.4 RS a 164 MHz pour champ_helio = 4 RS. ; tpmap = fltarr(npew, npew) ; usage limite a ce calcul. ; tpmap (rg(champ), sg(champ)) = 1 ; ; Remarque 1: le "placage" de la carte "champ", obtenue a par- ; ; tir des grilles xn et yn de l'image finale, sur ; ; la grille des coordonnees interferometriques rg ; ; et sg, n'est pas correct car ces grilles sont ; ; differentes. ; ; Remarque 2: rg et sg sont a-priori non-entiers. IDL accepte ; ; des indices fractionnaires pour tpmap (a condi- ; ; tion qu'ils soient dans [0, npew]) et les tron- ; ; que. ; ; Remarque 3: il peut arriver, a haute frequence notemment, que ; ; le choix de champ_helio = 4Rs par l'operateur ; ; fasse que les valeurs pour rg et sg debordent de ; ; [0, npew]. Leur pas de variation est alors > 1 ; ; et le remplissage de tmap peut etre lacunaire, ; ; d'ou la cuisine suivante pour boucher les trous. ;; Remplissage continu de l'aire de non-aliasing (cf Rem 3) ; for i = 0, npew - 1 do begin ; boucle sur les lignes. ; tab = where((tpmap(*, i) eq 1), c) ; ; tab contient les numeros des pixels non aliases de la ; ; ligne i de l'image. ; if (c gt 0) then tpmap(tab(0):tab(c-1), i) = 1 ; ; on decrete que tous les points non aliases sont conti- ; ; gus. ; endfor ; for i = 0, npew - 1 do begin ; boucle sur les colonnes. ; tab = where(tpmap(i, *) eq 1, c) ; if (c gt 0) then tpmap (i, tab(0):tab(c-1)) = 1 ; endfor ; ; champ = where(tpmap eq 1) ; ; champ est le drapeau de permission de calcul de champ, pre- ; ; tendument relatif a rg et sg. C'est un lonarr (re-indicage ; ; de tpmap, mettant les lignes bout a bout). ; endif ; fin du calcul des limites du champ facon J-M (contestable). ; if (keyword_set(selfga_np)) then harm_N = harm_N * selfga_np ; dans le cas d'une auto-calibration, les coeff selfga ont deja ; ete calcules dans une 1ere etape.