Pro RH_FIT_LIN_2D, af, nx, dnx, ny, dny, $	; entree
		bf				; sortie

; Creation 12 avril 2000.

; But : Extrapoler (une fonction de visibilite) sur les axes. On utilise un fit
;	  lineaire :
;		fit = a  +  b * i  +  c * j
;	  contraint sur le domaine (nx-dnx : nx , ny-dny : ny+dny)
; L'origine des coord est prise sur les axes au dernier point connu du champ de
;   contrainte.
; Le fit est limite au 1er degre car un fit polynomial de degre plus eleve 
;   divergerait tres vite a l'exterieur du domaine 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)

; Modifications:
; 00 avr 18	definition d'une grille de ponderation des points de af. Les
;		  premieres simulation sont decevantes: l'extrapolation n'est
;		  pas tres bonne surtout quand il y a une courbure et l'appli-
;		  cation a la visibilite donne pour les amul (dans CALIBRATION)
;		  des resultats plus mauvais que sans extrapolation. 
;		Le probleme est d'obtenir une meilleure evaluation du develop-
;		  pement limite au 1er ordre. On tente d'y arriver en ponderant
;		  davantage les points au bord du domaine de contrainte, la ou
;		  on extrapole.

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

; Notations:
; af		tableau d'entree, de dimensions (np * np). En pratique af est
;		  nul au dela de la bordure du champ 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)
    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_LIN_2D : tableau non carre en entree, a verifier'
	stop
    endif

    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 est un fltarr(dnx+1).

; Calcul de la moyenne
    moy = total (abs(fc)) / (dnx + 1)

; Recherche des indices des points contraignant le fit
    frac = 0.2				; pas trop grand car on peut devoir 
					;   extra poler pres de zeros de af.
;   tb = where(abs(fc) gt frac*moy)	; contient absc des pts admissibles.
    tb = where(abs(fc) ge    0)		; pour admettre tous les points.
    x  = tb - dnx			; origine au dernier point connu.

; Definition de la fonction de ponderation
	; les pts 0 et -1 ont meme poids. le pt -2 est deux fois moins lourd 
	;   que le pt -1, le point -3 est 3 fois mois lourd etc.
    w = fltarr(dnx+1)			; comme fc
    w(dnx)     = 1			; 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, ' x =', x
	print, ' w =', w
	print, 'fc =', fc

; Calcul des sommes du calcul des moindres carres.
    xt  = x (tb)
    fct = fc(tb)
    wt  = w (tb)
;   S0  = n_elements(tb)
    S0  = total (wt           )
    Sx  = total (wt * xt      )
    Sx2 = total (wt * xt  * xt)

    F0  = total (wt * fct     )
    Fx  = total (wt * fct * xt)

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

    vec = [F0,  Fx]

    ab = CRAMER(mat, vec)
    A  = ab(0)
    B  = ab(1)

; Vecteur des abscisses du champ de calcul du fit
    vx = indgen(dnx + 1 + 3) - dnx	; on prolonge de 3 pts sur l'axe (mais
					;   CALIBRATION 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.
	; 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
    Fit = A  +  B * vx

; Restriction du fit au 3 pts du champ de calcul qui excedent le champ de 
	; contrainte
    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 = 1
    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
	endif
	if (ny eq 0) then  begin
	    plot,  absc, bf(np/2 : np-1, np/2), xrange=[28, 35], linestyle=1
	    oplot, absc, af(np/2 : np-1, np/2)
	endif
	window, 1	&	wset, 1
	plot, absc, bf(np/2 : np-1, np/2), linestyle=1
	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 un sous en-
	;   semble de ce tableau, ou l'amplitude des harmoniques est non nulle
	;   et suffisante. Ce sous-tableau n'est en general pas carre.
	; On rappelle que le domaie de contrainte est (nx-dnx:nx, ny-dnu:ny)
    fc = af(np/2+nx-dnx : np/2+nx,  np/2+ny-dny : np/2+ny)
	; fc est un fltarr(2*dnx+1, 2*dny+1) Le point milieu est (dnx, dny).

; Calcul d'une moyenne representative. Dans l'utilisation de fit_lin_2d 
	; - soit nx=1  et ny=17, 20 ou 23 (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. On prend donc comme moyenne la moyenne sur l'axe concer-
	;   ne (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( abs(fc(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( abs(fc(*, 0)) ) / (dnx +1)
;   endif
    frac = 0.2			; pas trop grand car on peut devoir extrapoler
				;  au voisinage d'un zero de af.
    tb = where(abs(fc) 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

    x  = x  - dnx	; indices ramenes au centre du champ du fit.
    y  = y  - dny

; Calcul des sommes du calcul des moindres carres.
    S0 = n_elements(tb)

    Sx = total (x(tb))
    Sy = total (y(tb))

    Sx2 = total (x(tb) * x(tb))
    Sxy = total (x(tb) * y(tb))
    Sy2 = total (y(tb) * y(tb))

    F0 = total (fc(tb)        )
    Fx = total (fc(tb) * x(tb))
    Fy = total (fc(tb) * y(tb))

    mat = [[S0,  Sx,  Sy ],  $
	   [Sx,  Sx2, Sxy],  $
	   [Sy,  Sxy, Sy2]]

    vec = [F0,  Fx,  Fy]

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

; Construction de deux matrices, mx et my, contenant les abscisses et ordon-
	;  nees des points ou le fit sera calcule. Bien qu'on n'ait besoin 
	;  d'extrapoler les visibilites 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) - nx	; Origine sur l'axe x.
	vy = indgen(dny + 1 + 3) - dny	; On prolonge de 3 pts sur l'axe y,
					;   (mais CALIBRATION n'en utilise que
					;   2 par prudence pour extrapoler
					;   en y). 
					; Origine sur axe y au dernier pt connu
    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 vecteurs des colonnes et lignes (rangees par 
	;  valeurs croissantes) ou sont situes les points du champ de calcul 
	;  du fit, comptees a partir du centre du champ de contrainte.
    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, 'RH_FIT_LIN_2D : verif mx, my etc.'
	print, 'n1 =', n1
	help, mx, my
	stop
    endif

; 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
    Fit = A  +  B*mx  +  C*my

; Restriction du fit a la portion du champ de calcul qui excede le champ de 
	; contrainte
    if (nx eq 2) then fit_extrapol = fit(*           , dny+1:dny+3)  ; extrap y
    if (ny eq 0) then fit_extrapol = fitd(dnx+1:dnx+3,      *     )  ; academ.

; 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 = fltarr(dnx + 2*2 + 1, dny + 4*2 + 1)
	bf_r = fltarr(dnx + 2*2 + 1, dny + 4*2 + 1)
	jmax = 9 + dny			; nombre de lignes (pair).
	for i=0, 4+dnx do begin		; boucle sur les lignes
	    for j=0, jmax/2-1 do begin	; ---------  la moitie des lignes
		af_r(i, j) = af_red(i, jmax-1-j)
		bf_r(i, j) = bf_red(i, jmax-1-j)
	    endfor
	endfor
	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 l'axe y
    if (icon eq 1) then begin
	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
	oplot, bf_red, linestyle=1
	window, 1	&	wset, 1
	plot, bf_red
	stop
    endif

a2: j_m = 0

    bf = float(bf)

    end			; fin de RH_FIT_LIN_2D