C
C   Copyright (c) 1997 Silvano Bonazzola
C   Copyright (c) 1998 Jerome Novak
C
C    This file is part of LORENE.
C
C    LORENE is free software; you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation; either version 2 of the License, or
C    (at your option) any later version.
C
C    LORENE is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C    You should have received a copy of the GNU General Public License
C    along with LORENE; if not, write to the Free Software
C    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
C
C


	SUBROUTINE INSMTS(NDEG,JOB,ERR0,ERR1,DEN0,DEN1)
C
C		ROUTINE POUR L'INTERPOLATION ASTUCIEUSE (SMOOTH). 
C		(MINIMIZATION DE LA SOMME DU CARRE DE LA DERIVEE 2EME)
C
C		ARGUMENTS DE LA ROUTINE:
C
C		NDEG   = TABLEAU CONTENENT EN IMPUT LE NOMBRE DES POINTS
C			 DE L'INTERVALE OU  LA FONCTION EST DEFINIE (NDEG(1))
C			 ET LE NOMBRE DES POINTS OU LA FONCTION DOIT ETRE
C                        INTERPOLEE (NDEG(2)), EN OUTPUT NDEG(3) CONTIENT
C			 LE NOMBRE DES POINTS DE L'INTERVALLE UNION DE ERR0
C			 ET ERR1, DANS NDEG(4) IL-Y-A LE NOMBRE DES POINTS
C			 DE L'INTERVALLE ERR1
C
C		JOB= DRAPEAU: IL DOIT ETRE = 1 LORSQE ON APPELLE LA ROU-
C		          TINE LA PREMIERE FOIS OU LORSQUE ON CHANGE,NR1,MR1
C			 ,ERR0,ERR1 SI JOB.NE.1 ON EVITE DE RCALCULER LA 
C			  MATRICE DU SYSTEME AINSI QUE D'AUTRES TABLEAUX
C
C		ERR0	= TABLEAU IMPUT CONTENENT LES ABSISSES DE POINTS
C			  (EN ODRE CROISSANT) 
C			  OU LA FONCTION A INTERPOLER EST DEFINIE. EN OUTPUT
C			  ERR0 EST L'UNIONE ORDONNEE DE ERR0 (A L'ENTRE) ET
C			  DE ERR1. LORSQUE JOB.NE.1, ERR0 N'EST PAS UTILISEE
C			  ET PARCONSEQUENT IL RSTE INCHANGE
C
C		ERR1	= TABLEAU CONTENENT EN IMPUT LES ABSCISSES DES POINTS
C			  (EN ORDRE CROISSANT) OU IL FAUT INTERPOLER LA 
C		          FONCTION. EN OUTPUT IL-Y-A LA NOUVELLE GRILLE,
C			  C'EST A DIRE LA GRILLE FOORMEE PAR LES POINTS OU
C			  LA FONCTION DOIT ETRE INTERPOLEE ET LES POINTS
C			  DE ERR1 COINCIDENTS (A 1.E-6 PRES) AVEC LES
C			  DE ERR0. DANS CE DERTNIERS POINTS LA FONCTION N'EST
C			  PAS INTERPOLEE ET PREND SA VALEUR D'ORIGINE.
C			  LORSQUE JOB.NE.1, ERR0 N'EST PAS UTILISEE
C			  ET PARCONSEQUENT IL RSTE INCHANGE
C				
C		
C		DEN0	= TABLEAU CONTENENT EN IMPUT LES VALEURS DE LA
C			  FONCTION A INTERPOLER ET EN SORTIE LES VALEURES
C			  DE LA FONCTION INTERPOLEE SUR ERR0.(LA FONCTION
C			  A LE MEMES VALEURES SUR LES ANCIENS POINTS DE ERR0)
C
C		DEN1	= TABLEAU CONTENENT EN OUTPUT LES VALEURES INTERPOLEES
C			  DE LA FONCTION ET LES VALEURS D'ORIGONE DANS LES
C			  POINTS COMUNS DE ERR0 ET ERR1.
C			  
	IMPLICIT NONE
C
C
C $Id: insmts.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $
C $Log: insmts.f,v $
C Revision 1.2  2012/03/30 12:12:42  j_novak
C Cleaning of fortran files
C
C Revision 1.1.1.1  2001/11/20 15:19:30  e_gourgoulhon
C LORENE
C
c Revision 1.5  1998/09/30  12:53:36  novak
c le message d avertissement "tous les points coincident" a ete remplace
c par JOB = -1
c
c Revision 1.4  1998/09/22  15:19:46  novak
c N129 mis a 2049
c
c Revision 1.3  1998/09/22  15:15:19  novak
c rectification des arguments de LEQT1S (voir commentaires dans DGBSV)
c
c Revision 1.2  1997/05/23  11:40:35  hyc
c *** empty log message ***
c
C Revision 1.1  1997/03/19 19:24:57  eric
C Initial revision
C
C
C $Header: /cvsroot/Lorene/F77/Source/insmts.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $
C
C
	character*120 header
	data header/'$Header: /cvsroot/Lorene/F77/Source/insmts.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $'/

	INTEGER NDEG,NR,NR1,N129,LR,MR,MR1,J,JR,IPOS,J0,LR1,
     1	JR1,JR2,IPOS1,JM2,JM1,JP1,JP2,JOB,IPOS2,J1,JB1,JB2
C
	DOUBLE PRECISION
     1	ERR0,ERR1,DEN0,DEN1,BBB,AAA,AA,BB,CC,X12,X0,BQ1,DD,EE
C
	PARAMETER (N129=2049)
	DIMENSION ERR0(*),DEN0(*),DEN1(*),ERR1(*),
     1	BBB(N129,5),AAA(N129,3),AA(N129),BB(N129),CC(N129),IPOS(N129)
	DIMENSION DD(N129),EE(N129)
	DIMENSION IPOS1(N129),IPOS2(N129),NDEG(*)
C
	SAVE AAA,BBB,AA,BB,CC,X0,X12,BQ1,NR,MR,JR1,IPOS,IPOS1,IPOS2
C
	IF(JOB.NE.1) GO TO 993
C
C		CAS OU LA ROUTINE SERAIT APPELLEE POUR LA PREMIERE FOIS
C
C		REDIFINITION DE LA NOUVELLE VARIABLE r
C
	NR1=NDEG(1)
	MR1=NDEG(2)
	NR=NR1-1
	MR=MR1-1
	IF(ERR1(1).LT.ERR0(1).OR.ERR1(MR1).GT.ERR0(NR1)) THEN
	PRINT*,'ROUTINE D INTERPLOATION INSMTS: ON DOIT AVOIR'
	PRINT*,'ERR1(1).GE.ERR0(1) AINSI QUE ERR0(NR1).GE.ERR1(MR1)'
	PRINT*,'ERR0(1),ERR1(1)=',ERR0(1),ERR1(1)
	PRINT*,'ERR0(NR1),ERR1(MR1)=',ERR0(NR1),ERR1(MR1)
	STOP
	ENDIF
C
	IF(NR1+MR1.GT.N129) THEN 
	PRINT*,'ROUTINE INSMTS: DIMENSIONS LOCALES INSUFFISANTES:'
	PRINT*,'NR1=MR1,N129=',NR1+MR1,N129
	STOP
	ENDIF 
C
	X0=0
	X12=1.E-06*(ERR0(NR1)-ERR0(1))
C
C		 TOUS LES POINTS DE ERR0 ET DE ERR1 QUI DIFFERENT DE X12
C		 SONT CONSIDERES COINCIDENTS
C
	JR=0
	J0=1
C
	JB1=1
	JB2=1
	IF(ABS(ERR0(1)-ERR1(1)).LT.X12) JB1=-1
	IF(ABS(ERR0(NR1)-ERR1(MR1)).LT.X12) JB2=-1
	J=0
	DO LR=1,NR1+MR1
	IPOS(LR)=0
	ENDDO
C
	DO 999 LR=1,NR1-1
	JR=JR+1
	BBB(JR,1)=ERR0(LR)
	IPOS(JR)=LR
C
	DO J=J0,MR1
	IF(ERR1(J).GE.ERR0(LR)) THEN
	IF(ERR1(J).LT.ERR0(LR+1)) THEN	
	JR=JR+1
	BBB(JR,1)=ERR1(J)
	IPOS(JR)=0
	ELSE	
	J0=J
	GO TO 998
	ENDIF
	ENDIF
 777	CONTINUE
	ENDDO
  998	CONTINUE
  999	CONTINUE
C
	JR=JR+1
	BBB(JR,1)=ERR0(NR1)	
	IPOS(JR)=NR1
C
	DO LR=1,JR
	IPOS1(LR)=0
	ENDDO
C
C		ON VA IDENTIFIER LES POINTS DONT LES COORDONNEES 
C		DIFFERENT DE MOINS DE X12.
C
	J0=0
	J=0
	DO LR=1,JR-1
	LR1=LR+1
	IF(ABS(BBB(LR1,1)-BBB(LR,1)).GT.X12) THEN
	J0=J0+1
	ERR0(J0)=BBB(LR,1)
	IPOS(J0)=IPOS(LR)
	ELSE
	J=J+1
	IPOS1(J)=J0+1
C
	IPOS2(J)=MAX0(IPOS(LR1),IPOS(LR))	
	ENDIF
	ENDDO
C
	J0=J0+1
	ERR0(J0)=BBB(JR,1)
	IPOS(J0)=IPOS(JR)
C
	DO LR=1,J
	JP1=IPOS1(LR)
	IPOS(JP1)=-IPOS2(LR)
	ENDDO	
	DO LR=1,J0
	IPOS1(LR)=IABS(IPOS(LR))
	ENDDO
	IF(JB1.LT.0) IPOS(1)=-1
	IF(JB2.LT.0) IPOS(J0)=-J0
C
	JR1=J0
	JR=JR1-1
C
C		DANS ERR0 IL-Y-A LA GRILLE UNION DE ERR0 ET DE ERR1
C
C		 PREPARATION DES COEFFICIENTS PERMETTANT LE CALCUL DE
C		 LA DERIVE 2EME
C
	DO LR=2,JR
	BBB(LR,1)=2/(ERR0(LR+1)-ERR0(LR-1))
	ENDDO
C
	AA(1)=0
	BB(1)=0
	CC(1)=0
	AA(JR1)=0
	BB(JR1)=0
	CC(JR1)=0
C
	DO LR=2,JR
	AA(LR)=BBB(LR,1)/(ERR0(LR)-ERR0(LR-1))
	ENDDO
	DO LR=2,JR
	BB(LR)=-(1/(ERR0(LR)-ERR0(LR-1))+
     1	1/(ERR0(LR+1)-ERR0(LR)))*BBB(LR,1)
	ENDDO
	DO LR=2,JR
	CC(LR)=BBB(LR,1)/(ERR0(LR+1)-ERR0(LR))
	ENDDO
C
C		PREPARATION DE LA MATRICE ET DU 2EME MEMBRE DU SYSTEME
C
	DO J=1,5
	DO LR=1,JR1
	BBB(LR,J)=0
	ENDDO
	ENDDO
C
	DO LR=1,JR1
	DEN1(LR)=0
	ENDDO
C
	J0=0
	IF(IPOS1(2).GT.0) GO TO 997
	J0=1
	BBB(1,3)=BB(2)**2+AA(3)**2
	DEN1(1)=-AA(2)*BB(2)*DEN0(IPOS1(1))
	IF(IPOS1(3).GT.0) THEN
	DEN1(1)=DEN1(1)-(BB(2)*CC(2)+AA(3)*BB(3))*DEN0(IPOS1(3))
	ELSE
	BBB(1,4)=BB(2)*CC(2)+AA(3)*BB(3)
	ENDIF
	IF(IPOS1(4).GT.0) THEN
	DEN1(1)=DEN1(1)-AA(3)*CC(3)*DEN0(IPOS1(4))
	ELSE
	BBB(1,5)=AA(3)*CC(3)
	ENDIF	
  997	CONTINUE
C
	IF(JR.GE.5) THEN
	DO 996 JR=3,JR1-2
	IF(IPOS1(JR).GT.0) GO TO 995
	J0=J0+1
	JM2=JR-2
	JM1=JR-1
	JP1=JR+1
	JP2=JR+2
	BBB(J0,3)=CC(JM1)**2+BB(JR)**2+AA(JP1)**2	
	IF(IPOS1(JM2).GT.0) THEN
	DEN1(J0)=-AA(JM1)*CC(JM1)*DEN0(IPOS1(JM2))
	ELSE
	BBB(J0,1)=AA(JM1)*CC(JM1)
	ENDIF
	IF(IPOS1(JM1).GT.0) THEN
	DEN1(J0)=DEN1(J0)-(CC(JM1)*BB(JM1)+AA(JR)*BB(JR))*
     1	DEN0(IPOS1(JM1))
	ELSE
	BBB(J0,2)=CC(JM1)*BB(JM1)+AA(JR)*BB(JR)
	ENDIF
	IF(IPOS1(JP1).GT.0) THEN
	DEN1(J0)=DEN1(J0)-(BB(JR)*CC(JR)+AA(JP1)*BB(JP1))*
     1	DEN0(IPOS1(JP1))
	ELSE
	BBB(J0,4)=AA(JP1)*BB(JP1)+BB(JR)*CC(JR)
	ENDIF
	IF(IPOS1(JP2).GT.0) THEN
	DEN1(J0)=DEN1(J0)-AA(JP1)*CC(JP1)*DEN0(IPOS1(JP2))
	ELSE
	BBB(J0,5)=AA(JP1)*CC(JP1)
	ENDIF
C
  995	CONTINUE	 
  996	CONTINUE	
	ENDIF
	JR=JR1-1
	JM2=JR-2
	JM1=JR-1
	JP1=JR+1
	IF(IPOS1(JR).GT.0) GO TO 994
	J0=J0+1
	BBB(J0,3)=CC(JM1)**2+BB(JR)**2
	DEN1(J0)=-BB(JR)*CC(JR)*DEN0(IPOS1(JR1))
	IF(IPOS1(JM2).GT.0) THEN
	DEN1(J0)=DEN1(J0)-CC(JM1)*AA(JM1)*DEN0(IPOS1(JM2))
	ELSE
	BBB(J0,1)=CC(JM1)*AA(JM1)
	ENDIF
	IF(IPOS1(JM1).GT.0) THEN
	DEN1(J0)=DEN1(J0)-(CC(JM1)*BB(JM1)+AA(JR)*BB(JR))*DEN0(IPOS1(JM1))
	ELSE
	BBB(J0,2)=CC(JM1)*BB(JM1)+AA(JR)*BB(JR)
	ENDIF
  994	CONTINUE
C
C		LA MATRICE EST REDUITE DANS LE CAS OU IL-Y-AURAIT DES
C		" TROUS" DANS LES POINTS A INTERPOLER.
C
	DO LR=1,J0
	IF(BBB(LR,2).EQ.X0) THEN
	BBB(LR,2)=BBB(LR,1)
	BBB(LR,1)=0
	ENDIF
	IF(BBB(LR,4).EQ.X0) THEN
	BBB(LR,4)=BBB(LR,5)
	BBB(LR,5)=0
	ENDIF
	ENDDO 
C
	IF(J0.EQ.0) THEN
		JOB = -1
	RETURN
	ENDIF
C
C		CAS OU IL N'-Y-AURAIT QU'UN POINT OU LA FONCTION DOIT
C		INTERPOLEE
C
	IF(J0.EQ.1) THEN
	BQ1=1/BBB(1,3)
	DEN1(1)=DEN1(1)*BQ1
	ELSE
	CALL LEQT1S(BBB,J0,2,2,N129,DEN1,1,N129,0,AAA,J)
	ENDIF
C	
C		ON FORME LE TABLEAU OUTPUT DEN1
C
	J=0
	DO LR=1,JR1
	IF(IPOS1(LR).EQ.0) THEN
	J=J+1		
	ERR1(LR)=DEN1(J)
	ELSE
	ERR1(LR)=DEN0(IPOS1(LR))
	ENDIF
	ENDDO
C
C		ON FORME LE TABLEAU OUTPUT DEN0
C
	J=0
	J1=0
	DO LR=1,JR1
	IF(IPOS(LR).EQ.0) THEN
	J=J+1
	J1=J1+1
	DD(J1)=DEN1(J)
	ENDIF
	IF(IPOS(LR).LT.0) THEN
	J1=J1+1
	DD(J1)=DEN0(IPOS1(LR))
	ENDIF
	ENDDO
C
	DO LR=1,JR1
	DEN0(LR)=ERR1(LR)
	ENDDO
	DO LR=1,J1
	DEN1(LR)=DD(LR)
	ENDDO
C
C		ON FORME LE TABLEAU OUTPUT ERR1
C
	J=0
	DO LR=1,JR1
	IF(IPOS(LR).LE.0) THEN
	J=J+1
	ERR1(J)=ERR0(LR)
	ENDIF
	ENDDO
C	
C		LES NOMBRE DES POINTS DE LGRILLE ERR0 ET ERR1 SONT
C		MIS DANS LE TABLEAU NDEG
C
	NDEG(3)=JR1
	NDEG(4)=J1
C	
	RETURN
C
  993	CONTINUE
C
C		CAS JOB.NE.1
C
C---------------------------------------------------------------
C
	DO LR=1,JR1
	DEN1(LR)=0
	ENDDO
C
C		ON CALCULE LE 2EME MEMBRE DU SYSTEME
C
	J0=0
	IF(IPOS1(2).GT.0) GO TO 887
	J0=1
	DEN1(1)=-AA(2)*BB(2)*DEN0(IPOS1(1))
	IF(IPOS1(3).GT.0) THEN
	DEN1(1)=DEN1(1)-(BB(2)*CC(2)+AA(3)*BB(3))*DEN0(IPOS1(3))
	ENDIF
	IF(IPOS1(4).GT.0) THEN
	DEN1(1)=DEN1(1)-AA(3)*CC(3)*DEN0(IPOS1(4))
	ENDIF	
  887	CONTINUE
C
	JR=JR1-1

	IF(JR.GE.5) THEN
	DO 886 JR=3,JR1-2
	IF(IPOS1(JR).GT.0) GO TO 885
	J0=J0+1
	JM2=JR-2
	JM1=JR-1
	JP1=JR+1
	JP2=JR+2
	IF(IPOS1(JM2).GT.0) THEN
	DEN1(J0)=-AA(JM1)*CC(JM1)*DEN0(IPOS1(JM2))
	ENDIF
	IF(IPOS1(JM1).GT.0) THEN
	DEN1(J0)=DEN1(J0)-(CC(JM1)*BB(JM1)+AA(JR)*BB(JR))*
     1	DEN0(IPOS1(JM1))
	ENDIF
	IF(IPOS1(JP1).GT.0) THEN
	DEN1(J0)=DEN1(J0)-(BB(JR)*CC(JR)+AA(JP1)*BB(JP1))*
     1	DEN0(IPOS1(JP1))
	ENDIF
	IF(IPOS1(JP2).GT.0) THEN
	DEN1(J0)=DEN1(J0)-AA(JP1)*CC(JP1)*DEN0(IPOS1(JP2))
	ENDIF
C
  885	CONTINUE	 
  886	CONTINUE	
	ENDIF
	JR=JR1-1
	JM2=JR-2
	JM1=JR-1
	JP1=JR+1
	IF(IPOS1(JR).GT.0) GO TO 884
	J0=J0+1
	DEN1(J0)=-BB(JR)*CC(JR)*DEN0(IPOS1(JR1))
	IF(IPOS1(JM2).GT.0) THEN
	DEN1(J0)=DEN1(J0)-CC(JM1)*AA(JM1)*DEN0(IPOS1(JM2))
	ENDIF
	IF(IPOS1(JM1).GT.0) THEN
	DEN1(J0)=DEN1(J0)-(CC(JM1)*BB(JM1)+AA(JR)*BB(JR))*DEN0(IPOS1(JM1))
	ENDIF
  884	CONTINUE
C
	IF(J0.EQ.1) THEN
	DEN1(1)=DEN1(1)*BQ1
	ELSE
	CALL LEQT1S(BBB,J0,2,2,N129,DEN1,1,N129,2,AAA,J)
	ENDIF
C	
C		ON FORME LES TABLEAUX OUTPUT DEN0 ET DEN1
C
	J=0
	DO LR=1,JR1
	IF(IPOS1(LR).EQ.0) THEN
	J=J+1		
	DD(LR)=DEN1(J)
	ELSE
	DD(LR)=DEN0(IPOS1(LR))
	ENDIF
	ENDDO
C
	J=0
	J1=0
	DO LR=1,JR1
	IF(IPOS(LR).EQ.0) THEN
	J=J+1
	J1=J1+1
	EE(J1)=DEN1(J)
	ENDIF
	IF(IPOS(LR).LT.0) THEN
	J1=J1+1
	EE(J1)=DEN0(IPOS1(LR))
	ENDIF
	ENDDO
C
	DO LR=1,JR1
	DEN0(LR)=DD(LR)
	ENDDO
	DO LR=1,J1
	DEN1(LR)=EE(LR)
	ENDDO
	RETURN
C
  100	FORMAT(1X,10E10.3)
  101	FORMAT(1X,' ')
  200	FORMAT(1X,20I6)
	END	
C

