PRO RH_FIT_EXTRAPOL, i_r_im,  af_comp,  nx,  dnx,  ny,  dny,  $	; entree
		     bf						; sortie

; Creation 19 avril 2000, a partir de fit_lin_2d.pro.
; Modifications:
; 00 avr 26	entree de la visibilite complexe (pour pouvoir mieux calculer
;		  le drapeau de participation au fit d'apres le module) et du
;		  drapeau de choix de traitement de partie reelle ou imaginai-
;		  re i_r_im (argument supplementaire en entree).

; But : Extrapoler (une fonction de visibilite) sur les axes. On utilise un fit
;   quadratique complet :
;	  fit = a  +  b * i  +  c * j  + d * i^2  +  e * i * j  +  f * j^2
;   contraint sur le domaine (nx-dnx : nx , ny-dny : ny+dny), avec une pondera-
;   tion plus elevee au bord de l'intervalle de contrainte au dela duquel on 
;   extrapole.

; L'origine des coord est prise sur les axes au dernier point connu du champ de
;   contrainte.
; L'extrapolation est calculee le long de l'axe qui correspond a n* non nul, et
;   sur l'intervalle dn* correspondant a n*.
; La valeur dny=0 permet de reconnaitre le cas du fit 1D  :
;			Fit = A + B*i + c*i^2

; ATTENTION: la procedure est specifique au cas traite par autocalibration. Le
;	cas 2D est usuellement nx=2, dnx=4. Et le test if (nx eq 2) est utilise
;	plusieurs fois.

; Notations:
; af_comp	tableau d'entree complexe, de dimensions (np * np). En prati-
;		  que af_comp est nul au dela de la bordure du champ de con-
;		  trainte du fit.
; af		partie reelle ou imaginaire de af_comp (selon valeur de i_r_im)
; afm		module de af_comp (sert a choisir les points de contrainte du 
;		  fit.
; nx, ny	coordonnees du centre du champ de contrainte du fit.
; dnx, dny	2*dnx +1, 2*dny + 1 dimensions du champ de contrainte
; bf		tableau de sortie, de memes dimensions que af. Identique a af
;		  sauf sur la portion du champ de calcul qui excede le champ de
;		  contrainte di fit.

    taille = size   (af_comp)
    n_dim  = taille (0)
    npx    = taille (1)
    if(n_dim gt 1) then npy = taille(2)
    if (npx ne npy) then begin
	print, 'RH_FIT_EXTRAPOL : tableau non carre en entree, a verifier'
	stop
    endif

    afm = abs(af_comp)
    if(i_r_im eq 1) then  af = float    (af_comp)
    if(i_r_im eq 2) then  af = imaginary(af_comp)

    np = npx
    if(dny ne 0) then goto, a1		; vers le fit 2D

; Fit 1D. Le seul cas du fit purement 1D est celui de l'extrapolation de la vi-
	; sibilite au dela de 1600m de base EW pour etalonner E2: nx=32, dnx=3.
; Creation du sous-tableau utilise pour le fit.
    fc  = reform(af (np/2+nx-dnx : np/2+nx,  np/2))    ; fc   fltarr(dnx+1).
    fmc = reform(afm(np/2+nx-dnx : np/2+nx,  np/2))    ; fmc  fltarr(dnx+1).

; Calcul de la moyenne
    moy = total (fmc) / (dnx + 1)			; fmc est un module.

; Recherche des indices des points contraignant le fit
    frac = 0.3				; pas trop grand car on peut devoir 
					;   extrapoler pres de zeros de afm.
;   tb = where(fmc gt frac*moy)		; contient absc des pts admissibles.
    tb = where(fmc ge    0)		; admettre tous les points (essais).
    x  = tb - dnx			; origine au dernier point connu.

; Definition de la fonction de ponderation
	; le pts 0 (a l'extremite droite de l'intervalle de contrainte) pese 
	; 1.5 ou 2, le point -1 pese 1, le point -2 1/2, le point -3 1/3 etc.
    w = fltarr(dnx+1)			; comme fc
    w(dnx)     = 1.5			; par definition
    for i=1, dnx  do begin
	w(dnx - i) = 1. / i		; w est max au bout du tableau, au der-
					;   nier point connu.
    endfor
;   w(*) = 1

	print, 'RH_FIT_EXTRAPOL :'
	print, '   x =', x
	print, '   w =', w
	print, '  fc =', fc

; Calcul des sommes du calcul des moindres carres.
	; C'est x et et non x(tb) qui figure dans les sommes, car x est l'abs-
	; cisse des points retenus. Faute cherchee (dans le cas 2D) des heures
	; avec Anne le 21 avril 2000. Par contre il faut trier fc et w.
    fct = fc(tb)
    wt  = w (tb)
;   S0  = n_elements(tb)
    S0  = total (wt       )
    Sx  = total (wt * x  )
    Sx2 = total (wt * x^2)
    Sx3 = total (wt * x^3)
    Sx4 = total (wt * x^4)

    F0  = total (wt * fct       )
    Fx  = total (wt * fct * x  )
    Fx2 = total (wt * fct * x^2)

    mat = [[S0 ,  Sx ,  Sx2],  $
	   [Sx ,  Sx2,  Sx3],  $
	   [Sx2,  Sx3,  Sx4]]

    vec = [F0,  Fx,  Fx2]

    abc = CRAMER(mat, vec)
    A   = abc(0)
    B   = abc(1)
    C   = abc(2)

; Vecteur des abscisses du champ de calcul du fit
    vx = indgen(dnx + 1 + 3) - dnx	; on prolonge de 3 pts sur l'axe (mais
					;   autocalibration n'en utilise que 2
					;   par prudence pour extrapoler en y).
					; L'origine est au dernier point connu.

; Calcul de la fonction fit sur champ de contrainte + champ d'extrapolation, et
	; restriction du fit au 3 pts qui excedent le champ de contrainte
    Fit = A  +  B * vx  +  C * vx^2
    fit_extrapol = fit(dnx+1 : dnx+3)

; Construction de bf, prolongeant af hors du champ de contrainte. Gaffe lors
	; du prolongement d'un tableau 2D par un tableau 1D quand ce n'est pas
	; suivant une ligne.
    bf = af					; initialisation
    if (nx eq 2) then  begin			; extrapolation sur axe y.
	for iy = 1, 3   do begin
	    bf(np/2,  np/2 + ny + iy) = fit_extrapol(iy - 1)
	endfor
    endif
    if (ny eq 0) then  begin			; extrapolation sur axe x.
	for ix = 1, 3   do begin
	    bf(np/2 + nx + ix,  np/2) = fit_extrapol(ix - 1)
	endfor
    endif
	; car si le fit est 1D, les fonctions d'entree et de sortie sont 2D.

; Controle du resultat
    icon = 0
    if (icon eq 1) then begin
	window, 0	&	wset, 0
	absc = indgen(41)
	if (nx eq 0) then  begin
	    plot,  absc, af(np/2,  np/2 : np-1)
	    oplot, absc, bf(np/2,  np/2 : np-1), linestyle=1
	    xyouts, /normal, 0.5, 0.97, alignment=0.5, $
		'RH_FIT_EXTRAPOL : coupes NS de af ___, bf ...'
	endif
	if (ny eq 0) then  begin
	    plot, xrange=[25, 35], linestyle=1, $
;		yrange=[-0.25, -0.15], $
		   absc, bf(np/2 : np-1, np/2) 
	    oplot, absc, af(np/2 : np-1, np/2)
	    xyouts, /normal, 0.5, 0.97, alignment=0.5, $
		'RH_FIT_EXTRAPOL : coupes EW de af ___, bf ...'
	endif
	window, 1	&	wset, 1
	plot, absc, bf(np/2 : np-1, np/2), linestyle=1
	xyouts, /normal, 0.5, 0.97, alignment=0.5, $
		'RH_FIT_EXTRAPOL : coupes EW de bf ...'
	stop
    endif

    goto, a2				; sauter le fit 2D.

a1: j_m = 0		; fin du saut du fit 1D. 

; FIT 2D.
; Creation du sous-tableau utilise pour le fit. Le fit se fait sur les points
	;   d'un sous ensemble de ce tableau (sous ensemble qui n'est en gene-
	;   ral pas carre) ou l'amplitude des harmoniques est non nulle et suf-
	;   fisante.
	; On rappelle que le domaine de contrainte est (nx-dnx:nx, ny-dny:ny)
    fc  = af (np/2+nx-dnx : np/2+nx,  np/2+ny-dny : np/2+ny)	; p reel ou im
    fmc = afm(np/2+nx-dnx : np/2+nx,  np/2+ny-dny : np/2+ny)	; module
	; fc et fmc sont des fltarr(dnx+1, dny+1).

; Mise au point ; creation d'une fonction connue
    icon = 0
    if (icon eq 1) then begin
	A1 = 1.3  & B1 = 0.9  & C1 = 0.7  & D1 = 1.1  & E1 = 1.2  & F1 = -1.5
	vx = indgen(dnx + 1) - dnx/2	; Origine au dernier pt connu sur axe y
	vy = indgen(dny + 1) - dny
	    ; vx et vy sont les abscisses et ordonnees des points du champ de 
	    ;   contrainte, avec origine au dernier pt connu sur l'axe y.
	mx = vx # replicate (1, n_elements(vy))	; nl lignes   identiques.
	my = replicate (1, n_elements(vx)) # vy	; nc colonnes identiques.
;	print, 'Simulation de fc'
;	print, '   mx ='
;	print, mx
;	print, '   my ='
;	print, my
	fc = A1  +  B1*mx  +  C1*my  +  D1*mx^2  +  E1*mx * my  +  F1*my^2
;     Trace de fc simule
	window, 0	&	wset, 0
	shade_surf, fc
	stop
    endif

; fin de generation d'une fonction connue.

; Calcul d'une moyenne representative. Dans l'utilisation par autocalibration:
	; - soit nx=2  et ny=17, 19, 21   (extrapolation sur l'axe NS),
	; - soit ny=0  et nx=32           (----------------------- EW)
	; Sur aucun des axes EW ou NS la distribution des harmoniques n'est 
	;   lacunaire (alors qu'elle l'est en l'axe EW hors de cet axe). On 
	;   prend donc comme moyenne la moyenne sur l'axe concerne (x ou y).
;  Moyenne sur l'axe y. On rappelle que fc selectionne af de nx-dnx a nx (usuel
	; nx=2 et dnx=4 car les harm sont lacunaires en EW).
    if (nx eq 2) then begin
	moy = total(fmc(dnx/2, *)) / (dny +1)
    endif
;  Moyenne sur l'axe x. Le cas ne se presente pas pour le fit 2D dans l'auto-
	; calibration. Il est ici pour la forme.
;   if (ny eq 0) then begin
;	moy = total(fmc(*, 0)) / (dnx +1)
;   endif
    frac = 0.2			; pas trop grand car on peut devoir extrapoler
				;  au voisinage d'un zero de af.
    tb = where(fmc gt frac*moy)

; Recherche des indices de ligne (y) et de colonne (x) correspondant aux
	; indices contenus dans tb.
    n_c = dnx + 1		; nbre de colonnes du champ du fit.
    y   = tb / n_c
    x   = tb - y * n_c

; Choix de l'origine des coordoonnees (x et y) au dernier point connu sur l'axe
;  Cas ou on extrapole 2D sur l'axe x (non utilise dams autocalibration)
    if ((ny eq 2) and (dny eq 4)) then begin
	; bla-bla
    endif
;  Cas ou on extrapole sue l'axe y (seul cas 2D utilise dans autocalibration)
    if ((nx eq 2) and (dnx eq 4)) then begin
	x  = x  - dnx/2
	y  = y  - dny
    endif
    icon = 0
    if (icon eq 1) then begin
	print, 'tb =', tb
	print, 'x ='
	print, x
	print, 'y ='
	print, y
	stop
    endif

; Definition de la fonction de ponderation (adaptee au cas particulier d'ici).
	; le pts 0 (a l'extremite droite de l'intervalle de contrainte) pese 
	; 1.5 ou 2, le point -1 pese 1, le point -2 1/2, le point -3 1/3 etc.
    w = fltarr(dnx+1, dny+1)			; comme fc
    for i=0, dnx  do begin		; boucle sur les colonnes.
	for j=0, dny  do begin		; -------------- lignes.
	    if ((i eq dnx/2) and (j eq dny)) then begin	
					; dernier point connu sur l'axe y.
		w(i, j) = 1.3		; par definition
	    endif else begin		; tous les autres points
		ad = sqrt ((dnx/2 -i)^2 + (dny - j)^2)	
					; ad est la distance au dernier point 
					;   connu sur l'axe y.
		w(i, j) = 1 / ad
	    endelse
	endfor 
    endfor
;   w(*) = 1				; mise au point seulement.
    icon = 0
    if (icon eq 1) then begin
	print, 'w ='
	print, w
	window, 0	&	wset, 0
	shade_surf, w
	stop
    endif

; Calcul des sommes du calcul des moindres carres.
	; C'est x et et non x(tb)et y(tb) qui figurent dans les sommes, car x 
	;   et y sont les coordonnees des points retenus. Faute cherchee des 
	;   heures avec Anne le 21 avril 2000. Par contre il faut trier fc et w
    fct = fc(tb)
    wt  = w (tb)

    S0  = total (wt    )

    Sx  = total (wt * x)
    Sy  = total (wt * y)

    Sx2 = total (wt * x^2      )
    Sxy = total (wt * x   * y  )
    Sy2 = total (wt *       y^2)

    Sx3   = total (wt * x^3      )
    Sx2y1 = total (wt * x^2 * y  )
    Sx1y2 = total (wt * x   * y^2)
    Sy3   = total (wt *       y^3)

    Sx4   = total (wt * x^4       )
    Sx3y1 = total (wt * x^3 * y  )
    Sx2y2 = total (wt * x^2 * y^2)
    Sx1y3 = total (wt * x   * y^3)
    Sy4   = total (wt *       y^4)

    F0  = total (wt * fct            )
    Fx  = total (wt * fct * x        )
    Fy  = total (wt * fct       * y  )
    Fx2 = total (wt * fct * x^2      )
    Fxy = total (wt * fct * x   * y  )
    Fy2 = total (wt * fct       * y^2)

    mat = [[S0 ,  Sx   ,  Sy   ,  Sx2  ,  Sxy  ,  Sy2  ],  $
	   [Sx ,  Sx2  ,  Sxy  ,  Sx3  ,  Sx2y1,  Sx1y2],  $
	   [Sy ,  Sxy  ,  Sy2  ,  Sx2y1,  Sx1y2,  Sy3  ],  $
	   [Sx2,  Sx3  ,  Sx2y1,  Sx4  ,  Sx3y1,  Sx2y2],  $
	   [Sxy,  Sx2y1,  Sx1y2,  Sx3y1,  Sx2y2,  Sx1y3],  $
	   [Sy2,  Sx1y2,  Sy3  ,  Sx2y2,  Sx1y3,  Sy4  ]]  

    vec = [F0,  Fx,  Fy,  Fx2,  Fxy,  Fy2]

    res = CRAMER(mat, vec)
    A = res(0)
    B = res(1)
    C = res(2)
    D = res(3)
    E = res(4)
    F = res(5)

    icon = 0
    if (icon eq 1) then begin
	print, 'Valeurs entrees pour le test:'
	print, '   A =', A1, '    B =', B1,  '    C =', C1
	print, '   D =', D1, '    E =', E1,  '    F =', F1
	print, ' '
	print, 'Resultat du calcul des moindres carres :'
	print, '   A =', A, '    B =', B,  '    C =', C
	print, '   D =', D, '    E =', E,  '    F =', F
	stop
    endif
    
; Construction de mx et my, matrices des abscisses et ordonnees des points ou 
	; le fit sera calcule. Bien qu'on n'ait besoin d'extrapoler les visibi-
	; lites que sur les axes, on definit ici des matrices, dont une des 
	; dimensions est superieure a celle du champ de contrainte du fit.
;  Cas de l'extrapolation sur l'axe y (seul cas pratique pour le fit 2D).
    if (nx eq 2) then begin
	vx = indgen(dnx + 1    ) - dnx/2
	vy = indgen(dny + 1 + 3) - dny	
		; Origine sur l'axe y au dernier point connu sur l'axe y. 
		; On prolonge de 3 points sur l'axe y, (mais autocalibration 
		;   n'en utilise que 2 par prudence pour extrapoler en y). 
    endif
;  Cas de l'extrapolation sur l'axe x, non rencontre pour le fit 2D. Il figure
	; ici sous sa forme reellement 2D.
;   if (ny eq 0) then begin
;	vx = findgen(dnx + 1 + 3) - dnx	; on prolonge de 3 pts sur l'axe x,
					;   origine dernier pt connu sur axe x.
;	vy = findgen(dny + 1) - ny
;   endif
	; vx et vy sont les abscisses et ordonnees des points du champ de cal-
	;   cul, avec origine au dernier pt connu sur l'axe y.
    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) then begin
	print, 'verif mx, my'
	print, 'dnx =', dnx, '     dny =', dny
	help, mx, my
	print, 'mx ='
	print, mx
	print, ' '
	print, 'my ='
	print, my
	stop
    endif

; Calcul de la fonction fit sur champ de contrainte + champ d'extrapolation et
	; restriction du fit a la portion du champ de calcul qui excede le 
	; champ de contrainte
    Fit = A  +  B*mx  +  C*my  +  D * mx^2  +  E * mx * my  +  F * my^2
    if (nx eq 2) then fit_extrapol = fit(     *     , dny+1:dny+3)  ; extrap y
    if (ny eq 0) then fit_extrapol = fit(dnx+1:dnx+3,      *     )  ; academ.
    icon = 0
    if (icon eq 1) then begin
	fc_et = fit			; initialisation
	fc_et(*, *) = 0			; et mise a zero
	fc_et(0, 0) = fc		; fc etendue a tout le champ de calcul.
	window, 0	&	wset, 0
	shade_surf, reverse(fc_et, 2)		; presenter fc_et du bon cote.
	window, 1	&	wset, 1
	shade_surf, reverse(fit         , 2)
;	shade_surf, reverse(fit_extrapol, 2)
	stop
    endif

; Calcul des coord du coin inferieur gauche du domaine de prolongation dans le
	; champ de af
    if (nx eq 2) then begin		; cas extrapolation sur y.
	x_deb = np/2 - dnx/2
	y_deb = np/2 + ny +1
    endif
    if (ny eq 0) then begin		; extrapol sur x. Non rencontre en 2D.
	x_deb = np/2 + nx + 1
	y_deb = np/2 - dny/2
    endif

; Construction de bf, prolongeant af hors du champ de contrainte
    bf = af			; initialisation
    bf(x_deb, y_deb) = fit_extrapol 

; Controle
    icon = 0	; Trace 2D de la region extrapolee.
    if (icon eq 1) then begin
	af_red = (af(x_deb-2 : x_deb+dnx+2,  y_deb-dny-4 : y_deb+4))
	bf_red = (bf(x_deb-2 : x_deb+dnx+2,  y_deb-dny-4 : y_deb+4))
;     Retournement de la fonction (pour etre vue du bon cote par shade_surf)
	af_r = af_red	&	af_r(*, *) = 0		; initialisation
	bf_r = bf_red	&	bf_r(*, *) = 0
	taille = size(af_r)
	dim_x = taille(1)
	af_r = reverse(af_red, 2)
	bf_r = reverse(bf_red, 2)
;     Trace
	window, 0	&	wset, 0
	shade_surf, -congrid(af_r, 512, 512, cubic=-0.5)
	xyouts, /normal, 0.5, 0.97, alignment=0.5, 'fonction initiale' 
	window, 1	&	wset, 1
	shade_surf, -congrid(bf_r, 512, 512, cubic=-0.5)
	xyouts, /normal, 0.5, 0.97, alignment=0.5, 'fonction prolongee' 
	stop
    endif
    icon = 0	; Trace 1D sur les axe x et y
    if (icon eq 1) then begin
;     Trace sur l'axe x
	af_red = reform(af(np/2:np/2+nx+5, np/2))
	bf_red = reform(bf(np/2:np/2+nx+5, np/2))
	window, 0	&	wset, 0
	plot,  af_red, title='Coupe EW'
	oplot, bf_red - 0.02, linestyle=1		; 0.02 pour separer
	window, 1	&	wset, 1
	plot, bf_red
	stop
;     Trace sur l'axe y
	af_red = reform(af(np/2, np/2:np/2+ny+5))
	bf_red = reform(bf(np/2, np/2:np/2+ny+5))
	window, 0	&	wset, 0
	plot,  af_red, title='Coupe NS'
	oplot, bf_red - 0.02, linestyle=1
	window, 1	&	wset, 1
	plot, bf_red
	stop
    endif

a2: j_m = 0

    bf = float(bf)

    end		; fin de RH_FIT_EXTRAPOL.

; Poubelle:
; Calcul de la fonction fit sur champ de contrainte + champ d'extrapolation.
	; L'extrapolation (des essais ont montre que le fit peut s'ecarter sen-
	;   siblement de la fin de la fonction a extrapoler quand il y a une 
	;   courbure). On rappelle que le seul cas est celui de E2.
	; On assujetit le fit a passer a mi-chemin entre les 2 derniers points
	;   connus (d'abscisses -1 et 0), c'est a dire prendre leur valeur 
	;   moyenne moy pour vx = dnx - 0.5  .
	; A est alors tel que:
	;	A  +  B * (-0.5) = moy    soit  A = moy + 0.5 * B
;   A = (fc(dnx - 1) + fc(dnx)) / 2  +  0.5 * B
; Calcul de la fonction fit sur champ de contrainte + champ d'extrapolation.
	; L'extrapolation doit etre la plus adaptee aux variations sur l'axe.
	;   Des essais ont montre que le fit peut s'ecarter sensiblement de la 
	;   fin de la fonction a extrapoler quand il y a une courbure. On rap-
	;   pelle que le cas d'interet ici est celui de l'extrapolation sur
	;   l'axe y, pour etalonner NS18 a 23.
	; On assujetit le fit a passer a mi-chemin entre les 2 derniers points
	;   connus sur l'axe (d'ordonnees -1 et 0), c'est a dire prendre leur 
	;   valeur moyenne moy pour vy = dnx - 0.5  .   A est alors tel que:
	;	A  +  B * 0  +  C * (-0.5) = moy    soit  A = moy + 0.5 * C
;    moy = (fc(dnx/2, dny-1) + fc(dnx/2, dny)) / 2
;    A   = (fc(dnx/2, dny-1) + fc(dnx/2, dny)) / 2  +  0.5 * C