c            mslib.f
c            *******
c--------------------------------
      subroutine SNAME(iu,name,ier)
      character name*6,nx*6
1     format(a6)
      ier=0
c
      do n=1,10000
      read(iu,1,END=10)nx
      if(nx.eq.name)return
      enddo
 10   ier=-1
      return
      end
c----------------------------
c                -----
      subroutine ZERO3(ia,ib,ic,id,milsig,iprej,nt,lec,i1,i2,ksubs,
     1                jecri,ar,br,cr,item,nrjm,nrej,y1,y2,ifin)
      integer*2 lec(3000)
      dimension nrej(90),z1(3000),z2(3000),cr(10)
      ifin=0
      ar=0.
      br=0.
      y1=0.
      y2=0.
      if(ia.eq.0)return
c
      ifin=1
      if(iprej.ne.0)then
        if(milsig.eq.0)milsig=3000
        call REJEC3(lec,i1,i2,ia,ib,ic,id,milsig,iprej,nt,
c            ------
     1             jecri,cr,item,nrjm,nrej,ifin)
        if(ifin.ne.0)return
        y1=fzero(i1,nt,cr)
        y2=fzero(i2,nt,cr)
c          -----
        goto30
      endif
c
      if(milsig.ne.0)then
        do i=i1,i2
        z1(i)=lec(i)
        enddo
      call parafit(z1,i1,i2,milsig,z2)
c          -------
      endif
c
      xab=0.5*(ia+ib)
      ipiv=(ib-ia)/100+1
      sn=0.
      sy=0.
        do10 i=ia,ib,ipiv
        if(lec(i).eq.32000)goto10
        sn=sn+1.
          if(milsig.eq.0)then
          PIV=LEC(I)
          sy=sy+PIV
          goto10
          endif
        sy=amax1(sy,z2(i))
10      continue
      if(sn.eq.0.)return
      if(milsig.ne.0)sn=1.
      yab=sy/sn
c
      if((ia.eq.ic).and.(ib.eq.id))then
      ar=0.
      br=yab
      y1=yab
      y2=yab
      ifin=0
      goto30
      endif
c
      xcd=0.5*(ic+id)
      ipiv=(id-ic)/100+1
      sn=0.
      sy=0.
        do20 i=ic,id,ipiv
        if(lec(i).eq.32000)goto20
        sn=sn+1.
          if(milsig.eq.0)then
          PIV=LEC(I)
          sy=sy+PIV
          goto20
          endif
        sy=amax1(sy,z2(i))
20      continue
      if(sn.eq.0.)return
      if(milsig.ne.0)sn=1.
      ifin=0
      ycd=sy/sn
      ar=(ycd-yab)/(xcd-xab)
      br=yab-ar*xab
      y1=ar*i1+br
      y2=ar*i2+br
c
30    continue
      if(jecri.ne.0)then
35    format(' xab,xcd,yab,ycd, ar,br ='6e14.6)
      write(6,35)xab,xcd,yab,ycd,ar,br
36    format(' ia,ib,ic,id,milsig,iprej,i1,i2,ksubs,item,nrjm,y1,y2,'
     1       'ifin')
37    format(11i6,2f8.1,i2)
      write(6,36)
      write(6,37)ia,ib,ic,id,milsig,iprej,i1,i2,ksubs,item,nrjm,y1,y2,
     1           ifin
38    format(' lec avant/apres(si substitution)  tous les 10')
39    format(20i6)
      write(6,38)
      write(6,39)(lec(i),i=i1,i2,10)
      endif
c
      if(ksubs.eq.0)return
        do40 i=i1,i2
        if(lec(i).eq.32000)goto40
          if(iprej.ne.0)then
          kpiv=lec(i)-fzero(i,nt,cr)
c                       -----
          lec(i)=32000
          if(iabs(kpiv).lt.32000)lec(i)=kpiv
          goto40
          endif
        kpiv=lec(i)-ar*FLOAT(i)-br
        lec(i)=32000
        if(iabs(kpiv).lt.32000)lec(i)=kpiv
40      continue
      if(jecri.ne.0)write(6,39)(lec(i),i=i1,i2,10)
      return
      end 
c------------------------------
      subroutine DPMCAR(X,Y,P,ND,NT,COEF)
c                ------
c Ce sous-programme [ Schneider 75 ] calcule par une 
c methode de moindres carres en double precision le 
c polynome de NT termes (10 max) Y=F(X) associe a ND mesures.
C
C P=poids
c
      DOUBLE PRECISION C,F,PIVOT
      DIMENSION X(ND),Y(ND),P(ND),COEF(10),F(10),C(11,11)
      IF(NT.GT.10)NT=10
      NTP=NT+1
C
      DO 2 I=1,NT
      DO 2 L=1,NTP
2     C(L,I)=0.
C
      DO10 L=1,ND
      F(1)=1
        DO7 I=2,NT
C                   CHANGEMENT D'ECHELLE X
7       F(I)=F(I-1)*X(L)/1000.
      DO10 M=1,NT
        DO9 K=M,NT
9       C(K,M)=C(K,M)+F(K)*F(M)*P(L)
10    C(NTP,M)=C(NTP,M)+F(M)*P(L)*Y(L)
C
      DO20 M=1,NTP
      DO20 K=M,NTP
20    C(M,K)=C(K,M)
C
      DO104 I=1,NT
      PIVOT=C(I,I)
        DO103 M=1,NT
103     C(M+1,I)=C(M+1,I)/PIVOT
      DO104 K=1,NT
      IF(I.EQ.K)GOTO104
      PIVOT=C(I,K)
        DO105 M=1,NT
105     C(M+1,K)=C(M+1,K)-C(M+1,I)*PIVOT
104   CONTINUE
C
      COEF(1)=C(NTP,1)
      DO106 I=2,NT
106   COEF(I)=C(NTP,I)/(1000.**(I-1))
      RETURN
      END
c----------------------
C                -------
      SUBROUTINE PARAFIT(Y,I1,I2,L,Z)
      DIMENSION Y(1),Z(1)
        IF(L.EQ.0)THEN
          DO I=I1,I2
          Z(I)=Y(I)
          ENDDO
        RETURN
        ENDIF
C
        IF(L.GE.10000)THEN
        X=0.
          DO I=I1,I2
          X=X+Y(I)
          ENDDO
        X=X/FLOAT(I2-I1+1)
          DO I=I1,I2
          Z(I)=X
          ENDDO
        RETURN
        ENDIF
C
      SX2=0.
      SX4=0.
        DO N=1,L
        N2=N*N
        SX2=SX2+N2
        SX4=SX4+N2*N2
        ENDDO
        SX2=2*SX2
        SX4=2*SX4
      D=(2*L+1)*SX4-SX2**2
      IA=I1+L
      IB=I2-L
      IF(IB.LT.IA)RETURN
C
      DO100 I=IA,IB
      SY=0.
      SYX=0.
      SYX2=0.
      IL1=I-L
      IL2=I+L
        DO IP=IL1,IL2
        DI=IP-I
        SY=SY+Y(IP)
        SYX=SYX+Y(IP)*DI
        SYX2=SYX2+Y(IP)*DI*DI
        ENDDO
      A=(SY*SX4-SYX2*SX2)/D
      Z(I)=A
C
      IF(I.EQ.IA.OR.I.EQ.IB)THEN
      B=SYX/SX2
      C=((2*L+1)*SYX2-SY*SX2)/D
      ENDIF
C
      IF(I.EQ.IA)THEN
        DO IP=I1,IA
        DI=IP-IA
        Z(IP)=A+DI*(B+DI*C)
        ENDDO
      ENDIF
C
      IF(I.EQ.IB)THEN
        DO IP=IB,I2
        DI=IP-IB
        Z(IP)=A+DI*(B+DI*C)
        ENDDO
      ENDIF
C
100   CONTINUE
      RETURN
      END
c----------------------
	subroutine IV(p,km,lb,xi,xv,jabs)
	dimension p(41)
c
c xv en unite canal   1<xv<41
c
c  si lb=100 ponderation par (dp/dk)**2
c
c jabs=0  priorite aux plus fortes pentes
c      1  absorption seulement,priorite intensites faibles
c      -1 emission seulement,priorite intensites fortes
c si pas de solution  xi=xv=32000
c
      if(lb.eq.100)then
      sdp2=0.
      sxdp2=0.
      spdp2=0.
        do k=2,km-1
        dp2=(p(k+1)-p(k-1))**2
        sdp2=sdp2+dp2
        sxdp2=sxdp2+k*dp2
        spdp2=spdp2+p(k)*dp2
        enddo
      xi=spdp2/sdp2
      xv=sxdp2/sdp2
      return
      endif
c
	pte=0.
        yy=1000000.
        zz=0.
	xv=32000.
	xi=32000.
C
        IF((LB.EQ.0).AND.(JABS.EQ.1))THEN
        PS=100000.
        DO K=1,KM
        PS=AMIN1(PS,P(K))
        IF(PS.EQ.P(K)) KC=K
        ENDDO
        XI=P(KC)
        XV=KC
        RETURN
        ENDIF
C     
	ksol=0
	k1=lb+2
	  do100 k=k1,km
	  ka=k-lb-1
	  kb=ka+lb
	  d1=p(kb)-p(ka)
	  d2=p(kb+1)-p(ka+1)
	  if((d1*d2).gt.0.) goto100
	  if(d1.eq.d2)goto100
	  if((jabs.eq.1).and.(d1.gt.d2))goto100
	  if((jabs.eq.-1).and.(d1.lt.d2))goto100
        y=p(ka)+p(kb)+p(ka+1)+p(kb+1)
c
      if(jabs.eq.-1)then
      zz=amax1(y,zz)
      if(zz.eq.y)ksol=k
      goto100
      endif
c
      if(jabs.eq.1)then
      yy=amin1(y,yy)
      if(yy.eq.y)ksol=k
      goto100
      endif
c
 	  pt=p(kb+1)-p(kb)+p(ka)-p(ka+1)
	  pt=abs(pt)
	  pte=amax1(pt,pte)
	  if(pt.eq.pte)ksol=k
100	  continue
c
	if(ksol.eq.0)return
	ka=ksol-lb-1
	kb=ka+lb
	d1=p(kb)-p(ka)
	d2=p(kb+1)-p(ka+1)
	del=d2-d1
	if(del.eq.0.)return
	x=-d1/del
	xv=ka+float(lb)/2.+x
	xi=(p(ka)+p(ka+lb))*(1.-x)+(p(ka+1)+p(ka+1+lb))*x
	xi=0.5*xi
	return
	end
c-----------------------------
C                  ------
	SUBROUTINE REJEC3(LEC,I1,I2,IA,IB,IC,ID,MILSIG,IPAS,NT,
     1                   JW,CR,ITEM,NRJ,NREJ,IFIN)
	INTEGER*2 LEC(3000),IBON(3000)
	DIMENSION XR(3000),YR(3000),PR(3000),SIGMA(20),CR(10),
     1            NREJ(90)
2       FORMAT(' PLUS DE 30 POINTS A REJETER')
3       FORMAT(' PLUS DE 20 ITERATIONS POUR LA REJECTION')
	IFIN=1
      DO N=1,10
      CR(N)=0.
      ENDDO
      DO I=1,3000
      PR(I)=1.
      ENDDO
        NS=0
        NSM=0
	  DO10 ITER=1,20
          ITEM=ITER
	  IF(ITER.NE.1)EM=SIG*MILSIG/1000.
          NSP=NSM
          NS=0
	    DO4 I=IA,ID,IPAS
            X=I
            Y=LEC(I)
	    IBON(I)=0
            IF(Y.EQ.32000.)GOTO4
	      IF(ITER.NE.1)THEN
	      E=Y-FZERO(I,NT,CR)
C                 -----
	      IF(ABS(E).GT.EM)GOTO4
	      ENDIF
	    IBON(I)=1
            NS=NS+1
            XR(NS)=X
            YR(NS)=Y
4	    CONTINUE
        if(NS.lt.10) return
        NSM=NS
        CALL DPMCAR(XR,YR,PR,NSM,NT,CR)
C            ------
5       IF(NSP.EQ.NSM)GOTO20
C
C CALCUL DE SIG, ET CHARGEMENT DE NREJ
	SIG=0.
          DO NRJ=1,30
          NREJ(NRJ)=0
          ENDDO
        NRJ=0
        S=0.
        NS=0
	  DO6 I=IA,ID,IPAS
	    IF(IBON(I).EQ.0)THEN
            NRJ=NRJ+1
              IF(NRJ.GT.90)THEN
                IF(JW.EQ.1)WRITE(6,2)
                RETURN
              ENDIF
            NREJ(NRJ)=I
            GOTO6
            ENDIF
          S=S+1.
          NS=NS+1
	  SIG=SIG+(YR(NS)-FZERO(I,NT,CR))**2
C                       -----
6	  CONTINUE
	SIG=SIG/S
        SIG=SQRT(SIG)
        SIGMA(ITER)=SIG
          IF(JW.EQ.1)WRITE(6,12)SIG
C
10	CONTINUE
12      FORMAT(' SIGMA REJECTION='10F8.1)
        WRITE(6,12)SIGMA
          IF(JW.EQ.1)WRITE(6,3)
	RETURN
20	IFIN=0
          IF(JW.EQ.1)THEN
30        FORMAT('  CR(10),SIG  /  I REJETES  '/11E11.3)
32        FORMAT(2X30I4)
          WRITE(6,30)CR,SIG
          WRITE(6,32)NREJ
          ENDIF      
	RETURN
	END
C-----------------------------
      function fzero(i,nt,c)
c              -----
      dimension c(10)
      fzero=c(1)
      if(nt.lt.2)return
      do n=2,nt
      fzero=fzero+c(n)*i**(n-1)
      enddo
      return
      end
c---------------------------------------
