PRO RH_FIT_QUADRIC, i_stop, i_sys_lin, $ ; entree af, n1, ident=ident, $ ; entree. tcon, tnul, n_tnul, $ ; sortie. con pour "contrainte" A, B, C, D, fit_a ; sortie ; D, A, B, C, fit_a ; sortie avant 10 mars 2000. ; Procedures utilisant RH_FIT_QUADRIC au 12 dec 01 : ; - RH_MALC_IM_1D pour extrapolation parabolique de la composante continue ; - RH_MALC_IM_2D -------------------------------------------------------- ; Procedures comportant un appel a RH_FIT_QUADRIC delaiss'e au au profit de ; RH_FIT_QUARTIC12 dec 01 : ; - FIT_TF_1D (AUTOCAL_SUB) pour modeliser. ; - FIT_TF_2D --------------------------- ; Creation : 17 au 20 sept 99. ; Modifications : ; 00 mar 9 complet'e a une dimension ; oct 13 les points au voisinage de l'origine sont exclus du fit pour pou- ; voir traiter le cas des bases large. RH_FIT_QUADRIC est donc ; restreint au calcul de la composante continue. La partie 1D ; n'a pas ete modifiee et ne sert plus du tout. ; nov 10 tcon et tnul sont produits en sortie car ils sont necessaires a ; RH_FIT_CUBIC pour le fit d'une phase apres le fit d'un module. ; 01 mar 22 abandon de le methode de Cramer (plante souvent) au profit de la ; methode SVD pour resoudre les systemes lineaires n*n. ; jun 5 choix de la methode de resolition des systemes lineaires avec ; i_sys_lin (1 pour CRAMER, 2 pour LU, 3 pour SVD), sauf dans le ; cas 1D ou on resoud directement le systeme de Cramer seul. ; But : determine un fit quadratique sans termes du 1er degre. Etait des- ; tine a fitter une visibilite autour de son origine. Insuffisant ; et a ete remplace par: ; - fit quartic pour les modules des modeles centraux de visibilite ; (contraints sur sur un domaine entourant l'origine), ; - fit_cubic pour les phases des modeles centraux et decales. ; - fit_quadric_2 qui inclut les termes lineaires et destine a ; l'ecretage (creation 5 oct 00) . ; Depuis l'introduction de fit_cubic, fit_quartic et interpol (dans ; soustrac_base) fit_quadric ne sert plus qu'au calcul de la compo- ; sante continue dans malc_im_2d. Comme cette evaluation ne marche ; pas dans le cas d'une base large, a partir du 13 oct 2000 on exclut ; les bas harmoniques du fit et du calcul du niveau minimum d'un ; harmonique pour qu'il contraigne le fit. ; A est la composante continue. A B C coeff des termes de degre 2. ; Determine par interpolation parabolique af(np/2, np2) a partir ; de af(np/2-n1 : np/2+n1, np/2-n1 : np/2+n1). ; On calcule: ; - un fit isotrope: F_i = A + B*(i^2 + j^2) ; - un fit anisotrope F_a = A + B*i^2 + C*i*j + D*j^2 ; - un fit anisotrope F_a = A*i^2 + B*j^2 + C*i*j + D (avant 10-3-00) ; Les fits sont limites a un carre de (2*n1 + 1) * (2*n1 + 1) pts centre ; sur le centre du champ. ; On ne tient compte que des points qui sont superieurs a une certaine ; fraction de la moyenne des pts entourant le point central (definition ; du tableau tcon). Ceci pour fitter de facon realiste la visibilite ; dans le cas ou il y a plus d'un centre actif sur le soleil (l'ampli- ; tude pouvant alors etre fortement anisotrope (modulation selon la di- ; rection joignant les centres. ; Le calcul isotrope est conserve seulement pour comparaison. ; Si le tableau en entree est a une dimension le fit l'est aussi. B et C ; sont alors bidon. ; le mot-cle ident ne sert qu'au cas a une dimension (pour eviter de modi- ; fier la liste des argument du cas a deux dimensions. ; Notations: ; af fonction a fitter ; n1 nbre de pts de part et d'autre du centre du champ sur lequel ; porte le fit. ; ident identificateur du cas EW traite (de 1 a 6). ; tcon adresses des pts contraignant le fit (sortie pour FIT_CUBIC). ; tnul ---------------- ne ---- pas ----------------------------- taille = size(af) n_dim = taille(0) if(n_dim eq 0) then begin print, 'RH_FIT_QUADRIC : le tableau entree a la dimension 0' stop endif if (n_dim eq 1) then goto, a1 if (n_dim eq 2) then goto, a2 if(n_dim gt 2) then begin print, 'RH_FIT_QUADRIC : le tableau entree a une dimension > 2' stop endif ; CAS A UNE DIMENSION a1: J_M = 0 ; Generation d'une fonction connue imp0 = 0 if (imp0 eq 1) then begin n = 128 ad = findgen(n) - n/2 ; distance au point central (n/2). r2 = ad^2 amul = 0.3 ; variation entre le centre et n1. A0 = 1.0 B0 = -amul / n1^2 af = A0 + B0 * r2 plot, af(n/2-n1 : n/2+n1), $ xstyle=1, title='fit_quadric : fonction d essai' print, 'Fonction connue : A0=', A0, ' B0 =', B0 stop endif ; conclusion: ; Fin de controle npx = taille(1) n2 = npx / 2 ; Creation du sous-tableau utilise pour le fit. Le fit se fait sur un sous en- ; semble de ce tableau, ou l'amplitude des harmoniques est suffisante. ; La fonction entree fc est souvent un module de TF, donc symetrique ; par rapport au centre du tableau. On n'utilise pas ici cette propri- ; ete. if (ident ne 5) then begin ; le domaine du fit est d'un seul fc = af(n2-n1:n2+n1) ; tenant autour de l'origine endif if (ident eq 5) then begin ; le domaine du fit a un trou autour fc = af(n2-6-n1:n2+6+n1) ; de l'origine endif ; fc est un fltarr(2*n1+1). le centre est fc(n1) sauf pour ident=5 ou ; le centre est en fc(n1+6). ; Calcul d'une moyenne de l'amplitude des harmoniques. if ((ident eq 1) or (ident eq 3)) then begin ; harm EW pairs. if (imp0 eq 0) then begin ; vestige mise au point fc_moy = total(fc(n1-4:n1+4)) / 4 ; moy sur -4 -4 2 4 endif if (imp0 eq 1) then begin ; vestige mise au point fc_moy = total(fc(n1-1:n1+1)) / 3 ; pour fonction connue endif endif if (ident eq 2) then begin ; harm EW impairs. fc_moy = total(fc(n1-5:n1+5)) / 6 ; moy sur -5 -3 -1 1 3 5 endif if (ident eq 4) then begin ; harm NS avec NS0_ew. fc_moy = total(fc(n1-4:n1+4)) / 8 ; moy sur -4 a 4, 0 exclus endif if (ident eq 5) then begin ; harm NS avec NS0_ew. fc_moy = total(fc(n1+6-9:n1+6+9)) / 6 ; moy sur 7 a 9 endif if (ident eq 6) then begin ; harm NS avec NS0_ew. fc_moy = total(fc(n1-4:n1+4)) / 8 ; moy sur -4 a 4, 0 exclus endif frac = 0.3 ; fraction min de la moyenne pour um harm acceptable. fc_min = frac * fc_moy tcon = where(fc gt fc_min) ; adresses des pts ou fc > fc_min. tnul = where(abs(fc) lt fc_min, n_tnul) ; sortie pour RH_FIT_CUBIC. x = tcon - n1 ; abscisses des points fc > fc_min. n_x = n_elements(x) n2 = npx / 2 ; Controle de x icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_FIT_QUADRIC : controle de x' plot, fc, xstyle=1 print, ' tcon =', tcon print, ' x =', x stop endif ; conclusion: OK ; Fin de controle ; Calcul de sommes du calcul des moindres carres x2 = x^2 ; rappel: x est le vecteur des abscisses ou fc > fc_min. x4 = x^4 S0 = n_x S2 = total (x2) S4 = total (x4) F0 = total (fc(tcon)) ; total(fc) = total (fc(x)) F2 = total (fc(tcon) * x2) ; Calcul des determinants et des parametres de la parabole. Det = S4 * S0 - S2 * S2 DetA = S4 * F0 - S2 * F2 DetB = F2 * S0 - F0 * S2 AI = DetA / Det BI = DetB / Det c_c_p = AI ; Controle des sommes icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'S0 =', S0, ' S2 =', S2, ' S4 =', S4 print, 'F0 =', F0, ' F2 =', F2 print, 'Ai =', Ai, ' BI =', BI stop endif ; Construction du vecteur vx, de dimensions (2*n1+1), comme le tableau fc ; defini plus haut. Le modele est calcule sur vx, de dimension bien ; determinee, alors que le calcul du fit se fait sur le sous-ensemble ; de ce domaine ou les harmoniques ont des amplitudes suffisantes. vx = indgen(2*n1+1) - n1 ; Calcul du fit sur tout le domaine du fit (pts non fittables inclus). Fit = AI + BI * vx^2 dif = fc(tcon) - Fit(tcon) res = sqrt ( total(dif^2) / total(fc^2) ) ; residu relatif ; Controle du fit et visualisation. icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin ; Visualisation de la partie de af a fitter et du fit sur le meme domaine. window, 0 & wset, 0 plot, fc, xstyle=1, title='fonction en entree (domaine du fit)' window, 1 & wset, 1 plot, Fit, xstyle=1, title='fonction fit (domaine du fit)' if (imp0 eq 1) then begin print, 'Comparaison des valeurs vraies et fittees de A et B' print, ' A0 =', A0, ' B0 =', B0 endif print, 'Fit :' ch_format = "(' A =', f8.3, ' B =', f8.5)" print, format=ch_format, AI , BI print, ' residu relatif :', res stop endif ; Fin de controle ; Definition des parametres de sortie (notations commumes avec le cas 2D). A = AI B = BI C = 0 D = 0 fit_a = Fit goto, a3 ; FIN DU CAS A UNE DIMENSION ; CAS A DEUX DIMENSIONS a2: J_M = 0 ; Generation d'une fonction connue imp0 = 0 if (imp0 eq 1) then begin n = 128 ad = shift(dist(n), n/2, n/2) ; "distance", 0 en (n/2,n/2), n/2 au bord. r2 = ad^2 amul = 0.30 ; variation entre le centre et n1. A0 = 1.0 B0 = -amul / n1^2 af = A0 + B0 * r2 ; af(*, *) = 1 shade_surf, congrid(af(n/2-n1:n/2+n1, n/2-n1:n/2+n1), 512, 512, $ cubic=-0.5) window, 1 & wset, 1 tvscl, congrid(af(n/2-n1:n/2+n1, n/2-n1:n/2+n1), 512, 512, $ cubic=-0.5) stop endif ; conclusion: ; Fin de controle dim = size(af) npx = dim(1) npy = dim(2) if(npx ne npy) then begin print, 'RH_FIT_QUADRIC : le tableau d entree n est pas carre' print, ' Dimensions :', nx, ny stop endif n2 = npx / 2 ; Creation du sous-tableau utilise pour le fit. Le fit se fait sur un sous en- ; semble de ce tableau, ou l'amplitude des harmoniques est suffisante. fc = af(n2-n1:n2+n1, n2-n1:n2+n1) ; fc est un fltarr(2*n1+1, 2*n1+1). le centre est fc(n1, n1). ; Mise a zero de toutes les bases inferieures a 150m (13 oct 2000) fc(n1-3:n1+3, n1-3:n1+3) = 0 ; 49 elements mis a zero. ; Calcul d'un niveau minimum pour participer a la contrainte. fc_min = 0.2 * total(fc) / (n_elements(fc) - 49) ; car 49 elements a zero. ; fc_min = 0.0 * total(fc) / (n_elements(fc) - 49) ; pour essais. ; Calcul d'avant le 13 oct 2000. Calcul d'une moyenne. On somme a part les ; harm de base 50m et de base 100m. ; S_E0 = fc(n1+1, n1) + fc(n1-1, n1) ; S_50m = S_E0 + fc(n1, n1-1) + fc(n1, n1+1) ; S_100m = $ ; fc(n1+2, n1 ) + fc(n1+2, n1+2) + fc(n1 , n1+2) + fc(n1-2, n1+2) + $ ; fc(n1-2, n1 ) + fc(n1-2, n1-2) + fc(n1 , n1-2) + fc(n1+2, n1-2) ; frac = 0.3 ; fraction min de la moyenne pour um harm acceptable. ; if (S_E0 eq 0) then begin ; harm EW de 50m de E0 nuls. ; fc_min = frac * (S_50m + S_100m) / 10 ; endif ; if (S_E0 gt 0) then begin ; harm EW de 50m de E0 non nuls. ; fc_min = frac * (S_50m + S_100m) / 12 ; endif ; Calcul d'avant le 10 dec 99. ; fc_min = 0.3 * (fc(n1+2, n1 ) + fc(n1+2, n1+1) + $ ; reli'e a la moyenne ; fc(n1 , n1+1) + fc(n1-2, n1+1) + $ ; des pts pres de ; fc(n1-2, n1 ) + fc(n1-2, n1-1) + $ ; l'origine. ; fc(n1 , n1-1) + fc(n1+2, n1-1) ) /8 tcon = where(fc gt fc_min) ; adresses points fc > fc_min tnul = where(abs(fc) lt fc_min, n_tnul) ; sortie utile pout RH_FIT_CUBIC. ; Recherche des indices de ligne (y) et de colonne (x) correspondant aux ; indices contenus dans tcon. nr = 2 * n1 + 1 ; nbre de lignes et colonnes du champ reduit. y = tcon / nr x = tcon - y * nr x = x - n1 ; indices ramenes au centre du champ reduit. y = y - n1 n_x = n_elements(x) n_y = n_elements(y) if (n_x ne n_y) then begin print, 'RH_FIT_QUADRIC: nbre d abcisses different du nbre d ordonnees' print, ' n_x, n_y :', n_x, n_y stop endif n2 = npx / 2 ; Controle de x et y icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'fit_quadric : controle de x et y' window, 0 & wset, 0 tvscl, congrid(fc, 512, 512, cubic) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'tableau de contrainte du fit' ; print, 'x =', x ; print, 'y =', y window, 1 & wset, 1 plot, x, y, psym=2 stop endif ; conclusion: OK ; Fin de controle ; Passage en double precision x = long(x) & y = long(y) & fc = double (fc) ; Calcul de sommes du calcul des moindres carres solution isotrope r2 = x^2 + y^2 ; carre de la distance au centre du champ. r4 = r2^2 S0 = n_x S2 = total (r2) S4 = total (r4) F0 = total (fc(tcon)) ; total(fc) = total (fc(tcon)) F2 = total (fc(tcon) * r2) ; Calcul des determinants et des parametres du paraboloide de revolution. Det = S4 * S0 - S2 * S2 DetA = S4 * F0 - S2 * F2 DetB = F2 * S0 - F0 * S2 AI = DetA / Det ; A isotrope BI = DetB / Det c_c_p = AI ; Calcul de sommes du calcul des moindres carres solution anisotrope Sx4 = total(x^4) Sx2y2 = total(x^2 * y^2) Sx3y1 = total(x^3 * y ) Sx2 = total(x^2) Sy4 = total(y^4) Sx1y3 = total(x * y^3) Sy2 = total(y^2) Sx1y1 = total(x * y) S0 = n_x Fx2 = total(fc(tcon) * x^2) Fy2 = total(fc(tcon) * y^2) Fxy = total(fc(tcon) * x * y) F0 = total(fc(tcon)) mat = [[S0 , Sx2 , Sx1y1, Sy2 ], $ [Sx2 , Sx4 , Sx3y1, Sx2y2], $ [Sx1y1, Sx3y1, Sx2y2, Sx1y3], $ [Sy2 , Sx2y2, Sx1y3, Sy4 ]] vec = [F0, Fx2, Fxy, Fy2] ; Restauration de la simple precision x = long(x) & y = long(y) & fc = double (fc) ; Precaution pour tenter de reduire l'apparition du diagnostic de plantage ; "LUDC: Singular matrix in routine ludcmp" (ca reduit le determinant ; a 1) ; det_mat = determ(mat) ; mat = mat / det_mat^0.25 ; vec = vec / det_mat^0.25 ; Precaution upprimee car ca peut faire planter CRAMER (Input array is ; singular), alors que ca passe sans la normalisation du determinant. ; Controle : verif que la matrice d'entree n'est pas singuliere. icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, $ 'RH_FIT_QUADRIC, avant appel CRAMER, LINEAR_LU ou LINEAR_SVD: mat =' print, mat print, ' determ(mat) =', determ(mat) stop ; Rem : La precaution ci-dessus evite le plantage dans LUDC (Singular ; matrix in routine ludcmp), mais ca continue a planter dans ; CRAMER (Input array is singular), bien que determ(mat) calcule ; soit non nul. endif ; Resolution if (i_sys_lin eq 1) then abcd = CRAMER (mat, vec, /zero) if (i_sys_lin eq 2) then abcd = RH_LINEAR_LU (mat, vec, i_stop) if (i_sys_lin eq 3) then abcd = RH_LINEAR_SVD(mat, vec, i_stop) ; /double (mat et vec sont deja en double) et /column inutiles pour ; SVD. ; Faute trouvee le 11 dec 01 : 2 pour SVD et 3 pour LU ; meme faute dans RH_FIT_QUARTIC (mais pas dans RH_FIT_CUBIC). A = float(abcd(0)) ; avoir le resultat en simple precision. B = float(abcd(1)) C = float(abcd(2)) D = float(abcd(3)) icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_FIT_QUADRIC, apres appel resolution : mat =' print, mat print, ' determ(mat) =', determ(mat) print, ' Fit = A + B*x^2 + C*x*y + D*y^2 print, ' A, B, C, D :', A, B, C, D stop endif ; Sauvegarde du 10 mars 2000 (ordonner selon les puissances croissantes). ;;; Calcul de sommes du calcul des moindres carres solution anisotrope ; Sx4 = total(x^4) ; Sx2y2 = total(x^2 * y^2) ; Sx3y1 = total(x^3 * y ) ; Sx2 = total(x^2) ; Sy4 = total(y^4) ; Sx1y3 = total(x * y^3) ; Sy2 = total(y^2) ; Sx1y1 = total(x * y) ; S0 = n_x ; Fx2 = total(fc(tcon) * x^2) ; Fy2 = total(fc(tcon) * y^2) ; Fxy = total(fc(tcon) * x * y) ; F0 = total(fc(tcon)) ; mat = [[Sx4 , Sx2y2, Sx3y1, Sx2 ], $ ; [Sx2y2, Sy4 , Sx1y3, Sy2 ], $ ; [Sx3y1, Sx1y3, Sx2y2, Sx1y1], $ ; [Sx2 , Sy2 , Sx1y1, S0 ]] ; vec = [Fx2, Fy2, Fxy, F0] ; abcd = CRAMER(mat, vec) ; A = abcd(0) ; B = abcd(1) ; C = abcd(2) ; D = abcd(3) ; Construction de deux matrices, mx et my, de dimensions (2*n1+1, 2*n1+1), ; comme le tableau fc defini plus haut. Le modele est calcule sur mx ; et my, de dimension bien determinees, alors que le calcul du fit aux ; moindres carres ne se fait que sur le sous-ensemble de ce domaine ; ou les harmoniques ont des amplitudes suffisantes. ; Avant le 10 dec 99 mx et my etaient definies par rapport a ce sous- ; domaine, n'etaient pas necessairement de dimension 2*n1+1 ni meme ; carrees, alors qu'on supposait que c'etait le cas dans autocalibra- ; tion.pro pour l'immersion dans un tableau (np, np) ce qui produiait ; des erreurs. ; ancien calcul avant le 13 dec 99; mx et my n'avaient pas la m dim que fc. ; nc = max(x) - min(x) + 1 ; nbre de colonnes ; nl = max(y) - min(y) + 1 ; nbre de lignes ; vx = findgen(nc) + min(x) ; vy = findgen(nl) + min(y) vx = findgen(2*n1+1) - n1 vy = findgen(2*n1+1) - n1 ; vx et vy sont les vecteurs des colonnes et lignes (rangees par ; valeurs croissantes) ou sont situes les points du champ reduit, ; comptees a partir du centre du champ. mx = vx # replicate (1, n_elements(vy)) ; nl lignes identiques. my = replicate (1, n_elements(vx)) # vy ; nc colonnes identiques. ; Controle icon = 0 if ((icon eq 1) and (i_stop eq 1)) then begin print, 'RH_FIT_QUADRIC : verif mx, my etc. ' print, 'n1 =', n1 help, mx, my stop endif ; Calcul de la fonction fit isotrope sur tout le champ (pts non fittes inclus). Fit_i = AI + BI * (mx^2 + my^2) dif_i = fc(tcon) - Fit_i(tcon) res_i = sqrt ( total(dif_i^2) / total(fc^2) ) ; residu relatif ; Calcul de la fonction fit anisotrope sur tout le champ. Fit_a = A + B*mx^2 + C*mx*my + D*my^2 dif_a = fc(tcon) - Fit_a(tcon) res_a = sqrt ( total(dif_a^2) / total(fc^2) ) c_c_p = A ; Controle du fit et visualisation. icon = 0 if (icon eq 1) then begin ; Visualisation de la partie de af a fitter et des fits. print, 'RH_FIT_QUADRIC' window, 0 & wset, 0 shade_surf, congrid(fc, 512, 512, cubic) xyouts, /normal, 0.5, 0.97, alignment=0.5, 'fonction originale' window, 1 & wset, 1 shade_surf, congrid(Fit_i, 512, 512, cubic=-0.5) xyouts, /normal, 0.5, 0.97, alignment=0.5, 'fit isotrope' window, 2 & wset, 2 shade_surf, congrid(Fit_a, 512, 512, cubic=-0.5) xyouts, /normal, 0.5, 0.97, alignment=0.5, 'fit anisotrope' ; print, 'RH_FIT_QUADRIC : max(fit anisotope) = ', max(Fit_a) if (imp0 eq 1) then begin print, 'Comparaison des valeurs vraies et fittees de A et B' print, ' A0 =', A0, ' B0 =', B0 endif print, ' Fit isotrope : A + B * (x^2 + y^2)' ch_format = "(' A =', f8.5, ' B =', f8.3)" print, format=ch_format, AI , BI print, ' residu relatif isotrope:', res_i print, ' ' print, ' Fit anisotrope : A + B*x^2 + C*x*y + D*y^2 ' ch_format = "(' A =', f8.5, ' B =', f8.5, ' C =', " + $ "f8.5, ' D =', f8.3)" print, format=ch_format, A, B, C, D stop endif ; conclusion: ; Fin de controle a3: J_M = 0 ; fin du saut de traitement a deux dimensions. end ; fin de RH_FIT_QUADRIC