;+ *********************************************************************** ; NAME: ; CALCPOS ; ; PURPOSE: ; Cette procedure calcule les coordonnees heliographiques d'un point ; a partir des coordonnes instrumentales ; ; CATEGORY: ; NRH general ; ; CALLING SEQUENCE: ; CALCPOS, freq, ih, mn, r1, s1, $ ; input ; x, y ; output ; ; INPUTS: ; FREQ Frequence d'observation em MHz ; IH Heure de l'image ; MN Minutes " ; R1 position instrumentale EW ; S1 Position instrumentale NS ; ; KEYWORD PARAMETERS: ; Non ; ; OUTPUTS: ; X position heliographique EW ; Y position heliographique NS ; ; COMMON BLOCKS: ; Non ; ; SIDE EFFECTS: ; Ce module contient les procedures intrernes suivantes: ; OMEG_BRUT calcul de l'angle solide du champ ; POS_HELIO calcul des positions ; HEURE_TU calcul du temps TU pour les sources ; ; MODIFICATION HISTORY: (A Kerdraon) ; Adapte de XHELIO (bonmartin@obspm.fr) ; 02 avr 3: mise au clair de commentaires sur new, nns, ins45 etc. dans ce ; fichier et dans INIT_MALAX_2D.PRO. ;-******************************************************************* ;_____________________________________________________________________________ FUNCTION OMEG_BRUT, freq, ih, mn ; But : calcul de l'angle solide du champ brut (1 periode dans chaque sens) 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 ah = !pi/12.* (ih + float(mn) / 60. - ahmer); soleil ou non solaire. cosh = cos (ah) sinh = sin (ah) ; Calcul de l'angle solide du champ (brut) al = c / freq ; longueur d'onde en metres. sinu = sinh * cosd cosu = sqrt(1 - sinu^2) Pew = al / (dew * 2 * cosu) ; periode ew fixee par la base de 100m sinu = cosh * cosd * sinl - sind * cosl cosu = sqrt(1 - sinu^2) Pns = al / (dns * cosu) omega = Pew * Pns return, omega end ; fin de OMEG_BRUT. ;______________________________________________________________________________ PRO POS_HELIO, freq, ah, r, s, $ ; entree x1, y1, x, y ; sortie: pos helio en min d'arc et RS. ; But : calculer les positions heliographiques a partir des positions interfe- ; rometriques r et s, lesquelles sont comprises en entree comme des ; 1/new et 1/nns des champs. ; corriger (eventuellement les erreurs de temps meridien (ou d'horloge) ; et de declinaison faites dans le programme d'observation temps reel. ; Il est prevu de corriger ces delouiseries historiques, mais les va- ; leurs de ces erreurs sont ici prises nulles. ; Remarque : les coordonnees heliographiques tiennent compte : ; . des derives en angle horaire et declinaison du soleil ; . de la rotation -p si rot_p=1 . Si rot_p=0, INIT_MALAX_2D pose ; sinp=0 et cosp=1 dans le common MALAX. ; ; Modifications : ; 02 avr 3: amelioration des commentaires sur dew, dns, new, nns dans ce ; fichier et dans INIT_MALAX_2D.PRO. ; Commentaires ; new et nns sont les nbres de points en EW et en NS sur l'image interfero ; metrique. Ce ne sont pas necessairement les nbres de points en puissance ; de 2 minimaux pour echantillonner correctement les images. ; A) Pour les images 2D (procedures en aval de RH_DPNEW) on a choisi de ; fixer dans RH_DPNEW new=nns=128 (not'es dans cette procedure npi, ; "nbre de pts sur l'image interferometrique"). Il faut donc, our uti- ; liser RH_DPNEW -> RH_DPATCHFITS_NRH -> RH_CALC_IMAGE_HELIO, remplir ; le common MALAX avec les valeurs de new et nns choisies. C'est ce ; qui est fait depuis le 10 sept 2001 en levant, dans INIT_MALAX_2D, ; les drapeaux ie0, ie1, ie2, ins45. => new et nns valent 128. ; B) Pour des images de l'ancien soft 2*1D, obtenues sans E0 ni NS45 et ; definies sur 64 points, il faut appeler IN_MALAX avec un autre ; appelant que INIT_MALAX_2D. ; Principe du calcul ; On calcule (en utilisant les parametres contenus dans le common MALAX) : ; . les angles horaires hew et hns, ; . lambda ; . les corrections aux r et s mesures, dues a l'ionosphere spherique et a ; l'effet de coin (ce qui est prevu en amont dans CALCPOS mais qui n'est ; pas fait realise au 18 sept 01) ; . les coordonnees interferometriques reduites (sans dimensions et compri- ; ses entre -1/2 et +1/2) : ; r / new * lambda / dew ; s / nns * lambda / dns ; . les ecarts angulaire cosddh et ddelta et leur conversion en RS ; . les coordonnees heliographiques dans un repere se deduisant du repere ; (H, delta) par une rotation de l'angle p (dans ce sens). ; Rappel algebrique (18 sept 01) ; On appelle : ; V_u vecteur unitaire reperant un point sur le soleil ; V_0 ------------------------- le centre du champ ; Dew vecteur de base du reseau EW, de module dew ; Dns ------------------------- NS----------- dns ; new nombre de points en EW sur l'image interferometrique ; nns ------------------- NS ----------------------------- ; ; On pose : delta(V_u) = V_u - V_0 (1) ; Les images interferometriques (qui sont simplement les TF (directes au sens ; d'IDL, avec un drapeau -1) des harmoniques mesures par le reseau) sont ; definies en new et nns points. Par hypothese ces points sont en nombre ; suffisant pour permettre d'echantillonner l'image et ils sont repartis sur ; une periode lambda/dew et lambda/dns sur chaque axe. ; Les relations qui permettent le passage des coordonnees interferometriques ; r et s, qui sont les rangs des points ou l'image est calculee, aux coor- ; donnees heliographiques sont : ; delta(V_u) . Dew = r * lambda / new (2 - 1) ; delta(V_u) . Dns = s * lambda / nns (2 - 2) ; Les equations (2) sont le systeme de depart du calcul litteral du 18 dec 90 ; qui inclut aussi les effets ionospheriques et qui sert de base a tous les ; programmes de calcul de position. ; Remarques : ; . les coordonnees heliographiques (en Rs) interviennent dans V_u. ; . new et dew interviennent par leur produit dans (2). De meme nns et dns. ; Il en resulte les points suivants : ; - l'introduction de E2 a necessite un changement des parametres ; utilises dans le calcul (new=32, dew=100m -> new=64, dew=100m), ; les points calcules sur l'image etant 2 fois plus serres). ; - l'introduction de E0 a ete transparente (les calculs de positions ; peuvant etre faits avec les anciennes valeurs (64 pts et 100m), ; les points etant 2 fois plus nombreux sur un champ 2 fois plus ; large. ; Le point nouveau avec E0 est que le champ lambda/dew etant dou- ; bl'e, le rang r peut etre un entier 2 fois plus grand qu'avant et ; peut decrire des sources loin du centre. ; - l'introduction des antennes impaires du NS a aussi ete transparente ; (nns=32, dns=108m -> nns=64, dns=54 m). ; - l'introduction de NS45 a necessite un changement (nns=64, dns=54m ; -> nns=128, dns=54m). ; . les nombres de points new et nns dans (2) ne correspondent pas necessai- ; rement a l'echantillonnage minimum (ou a la puissance de 2 immediate- ; ment superieure). Dans les faits cela a ete le cas avant la construc- ; tion de E0 et NS45 et apres l'utilisation de ces dernieres, mais il ; y a eu une periode intermediaire (depuis le debut du RH numerique vers ; 96 jusqu'a l'utilisation de NS45 vers jan 01) ou les images ont ete ; surchantillonnees a 128 points en EW et en NS. ; Cela a cr'e'e un certaine confusion et des erreurs. ; . l'image etant centree, r et s varient de -new/2 et -nns/2 a new/2-1 et ; nns/2-1. common MALAX, $ ; rempli par INIT_MALAX ; Voir projet de rearrangement du common MALAX dans INIT_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 sinhew = sin(ah - hew) coshew = cos(ah - hew) sinhns = sin(ah - hns) coshns = cos(ah - hns) dhc = 2 * !pi / 86400.* etmer ; erreur sur le temps meridien(radians) ddc = 2 * !pi / 360. * (edec / 60) ; ----la declinaison------- alambda = c / freq ; long d'onde en metres. ; Correction des coordonnes interferometriques ("canaux") r et s tenant compte ; des erreurs sur le temps meridien et la declinaison pendant l'obser- ; vation. Les canaux sont des 1/new et 1/nns des champs EW et NS. aew = new * dew / alambda ans = nns * dns / alambda dr = aew * (-cosd * coshew * dhc + (sind * sinhew -u * cosd) * ddc) ds = ans * ( cosd *sinhns * sinl * dhc + $ (sind * coshns *sinl + cosd * cosl) * ddc) r = r - dr s = s - ds ; print, 'corrections erreurs temps reel : dr =', dr, ' ds =', ds ; Calcul des ecarts angulaires cosddh et ddelta (en radians) entre le sursaut ; radio et le centre du disque. Correction de la derive du soleil en ; angle horaire et declinaison. ; Remarque: le calcul des positions heliographique ne fait intervenir ; que les quantites r/(new*dew) et s/(nns*dns). Les va- ; leurs utilisees ici pour new et nns sont lues dans le ; common MALAX et fixees en amont par INIT_MALAX (en general ; 128*64). Cela posera un probleme pour le calcul par ; MAILLE_3 des coord helio des 3 points definissant la maille ; elementaire d'une grille interferometrique de npew*npew ; (en general 128*128) points. r = r * c / freq / (dew * new) s = s * c / freq / (dns * nns) denom = sind * sinl * cos(hns - hew) + cosd * (cos(ah - hew) * cosl - $ u * sin(ah - hns) *sinl) cosddh = r * (sinl * sind * cos(ah - hns) + cosl * cosd) + $ s * (sin(ah - hew) * sind - u * cosd) ddelta = - r * sinl * sin(ah - hns) + s * cos(ah - hew) cosddh = cosddh / denom ; en radians. ddelta = ddelta / denom ; ---------- cosddh = cosddh + ghmer * ah * cosd ddelta = ddelta - gdel * ah ; Rotation d'angle -p ( p = (axe terrestre, axe solaire) ) et conversion en ; rayons solaires. ; Si rot_p=0 en entree dans RH_DPNEW, sinp et cosp ont et mis a 0 et ; 1 dans le common MALAX par INIT_MALAX_2D. ; x1 = cosddh * cosp + ddelta * sinp ; en radians y1 = -cosddh * sinp + ddelta * cosp ; ---------- x1 = x1 * 180 / !pi * 60 ; en min d'arc. y1 = y1 * 180 / !pi * 60 ; ------------ x = x1 / rsol ; en RS. y = y1 / rsol return end ; fin de POS_HELIO. ;______________________________________________________________________________ PRO HEURE_TU, isour, ian1, im, ij, ih, mn,ihtu, imtu ; ihtu et imtu sont utilises pour calculer l'etat de l'ionosphere: ; - identiques a ih et mn si isour = 1 (soleil). ; - differents de --------------------- 0 ( source). mois = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] if (isour eq 1) then begin ; soleil. ihtu = ih imtu = mn return endif ; Transformation du TS en TU: ; Calcul du TS de Greenwich a 0 h TU pour la date consideree ; On utilise l expression : ; TS Grenw. 0 h TU = 6h 41m 50.55sec + 8640184.81sec * t ; + 0.09sec * t**2 ; ou t est en siecles de 36 525 jours a partir de j2000.0 ; (Aoki, Guinot, Kaplan, Kinoshita, Mccarthy, Seidelman, 1982, ; Astron. & Astrophys. 105, 359) ian = rh_date100 ( ian1, /full) ; Calcul de t ta = (ian - 2000.) * 365. tm = 0. for i = 0, im - 2 do tm = tm + mois (i) t = ta + tm + ij - 1 ; en jours. ; Traitement des annees bissextiles. A VERIFIER (vu le 17 sept 01) : if (ian lt 1996) then t = t - 1 ; t<0 et un jour de plus if (ian eq 1996) and (im lt 3) then t = t - 1 if (ian lt 1992) then t = t - 1 if ((ian eq 1992) and (im lt 3)) then t = t - 1 if (ian lt 1988) then t = t - 1 if ((ian eq 1988) and (im lt 3)) then t = t - 1 if ian lt 1984 then t = t - 1 if ((ian eq 1984) and (im lt 3)) then t = t - 1 if (ian lt 1980) then t = t - 1 if ((ian eq 1980) and (im lt 3)) then t = t - 1 ; passage an 2000 if (ian gt 2000) then t = t + 1 if (ian eq 2000) and (im ge 3) then t = t + 1 if (ian gt 2004) then t = t + 1 if (ian eq 2004) and (im ge 3) then t = t + 1 if (ian gt 2008) then t = t + 1 if (ian eq 2008) and (im ge 3) then t = t + 1 t = t / 36525. ; en siecles. ; Calcul du TS de Greenwich a 0 h TU: ts0 = 6.664481 ; j ai retire 1min 58.36 sec pour que ca marche ts1 = 2400.051336 ts2 = 0.000026 tsG0tu = ts0 + ts1 * t + ts2 * t^2 ; TS de Green.a 0 h TU (heures) ; Cadrage entre 0 et 24 h : izit = fix (tsG0tu / 24) if (tsG0tu lt 0) then tsG0tu = tsG0tu + 24 * (-izit + 1); jusqu'en 2000. ; en heures. if (tsG0tu gt 0) then tsG0tu = tsG0tu + 24 * (-izit) ; apres 2000. ; en heures. ; Calcul du TS a Nancay a 0 h TU : along = - 523.3 ;longit. Nancay (sec, >0 a l'W). tsNan0tu = 3600 * tsG0tu - along ;TS de Nancay a 0 h TU (en sec). ; Calcul du TU en fonction du TS du lieu d'observation et cadrage 0 - 24 h. ahts = 3600 * ih + 60 * mn ; heure TS locale (sec). ahtu = 0.9972696 * (ahts - tsNan0tu) ; 0.9972696=365.25/366.25 if (ahtu lt 0) then ahtu = ahtu + 24 * 3600 if (ahtu gt 0) then ahtu = ahtu - 24 * 3600 ihtu = fix (ahtu / 3600) imtu = fix ( (ahtu - 3600 * ihtu) / 60) return end ; fin de HEURE_TU. ;______________________________________________________________________________ PRO CALCPOS, freq, ih, mn, r1, s1, $ ; entree x, y ; sortie ; But : calcul de positions, a appeler a chaque changement d'heure et de ; frequence. ; Commentaire detaille : voir plus haut dans POS_HELIO. ; Notations : ; r1 s1 coordonnees interferometriques a transformer. ; print, 'entree CALCPOS' ; print, 'freq', freq ; print, 'heure', ih, mn ; print, 'posit', r1, s1 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 ; Mise au point imp = 0 ; 0 sauter, 1 mise au point if (imp eq 1) then begin print, $ 'CALCPOS: new=', new, ' dew=', dew, ' nns=', nns, $ ' dns=', dns stop endif HEURE_TU, isour, ian, im, ij, ih, mn, ihtu, imtu ; ihtu et imtu sont utilises pour calculer l'etat de l'ionosphere: ; - identiques a ih et mn si isour = 1 (soleil). ; - differents de --------------------- 0 ( source). ; Calcul de angles horaires. ahsol = (ihtu + float(imtu) / 60. - ahmersol) ; soleil (heures). ahsol = !pi / 12.* ahsol ; ------ radians coshsol = cos(ahsol) sinhsol = sin(ahsol) ah = !pi/12.* (ih + float(mn) / 60. - ahmer) ;soleil ou non solaire. ; ah angle horaire calcule sans tenir compte des derives du soleil en ; ascension droite et declinaison (simplement 1 tour en 24 h TU. cosh = cos (ah) sinh = sin (ah) ; ici truc ionosphere a ajouter (appel cor_iono) r = r1 s = s1 POS_HELIO, freq, ah, r, s, $ ; entree x1, y1, x, y ; sortie: pos helio en ' d'arc et RS. ; variable dim pas initialisee si pas de ionosphere... ; dirs = dim / rsol ; distz dist zenithale au pt iono I (deg). ; dim deflection ionospherique (min d'arc). return end ; fin de CALCPOS ;_____________________________________________________________________________