        OPTIONS(DREAL)
c,OPTIMISE)
        PROGRAM EQWREC
C                       VERSION 3.2.2001
C       COMPILATION:  FTN77 EQWREC2 /LINK77
C
C       A CODE FOR:
C       -EQUIVALENT WIDTH MEASUREMENT
C       -RADIAL VELOCITY MEASUREMENT
C       -RECTIFICATION OF SPECTRA (SCALING TO LOCAL CONTINUUM=1)
C       -REMOVING A TREND (SUBSTRACTION OF A HAND-DEFINED FUNCTION)
C       -HAND-FITTING OF A FUNCTION
C       -ADDING SOME EQUIDISTANT X-VALUES TO SINGLE COLUMN DATA
C
c       authors: Budaj J., Komzik R.
c       contact: budaj@ta3.sk
c       published: http://www.ta3.sk/~budaj/software
c       special requirements: DBOS
c       non original routins:  intep (Press.et.al. 1992, Numerical recipes)
c       input: keyboard + file
c       output:  fit+file
c
        PARAMETER (NBOD=20000,NBGR=200)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION XP(NBGR),YP(NBGR)
        common /a/x(nbod),y(nbod)
        character*14 in,out
        character KEY*1
        character*30 curdir@,ipa
        logical ex

        open(4,file='FIT',status='unknown')
        call clear_screen@
        call hide_cursor@
C        call coup@(' W E L C O M E   T O ',116,29,20)
C        call coup@(' R E C T I F 7 7 ',244,31,21)
        ipa=curdir@()
        call coup@(' Path : => ',30,0,2)
        call read_edited_line@(ipa,11,2,113,ic)
        call attach@(ipa,ier)
        call doserr@(ier)
10      call coup@(' Input file name : => ',30,0,4)
16      call select_file@('*.*',in,*16)
        call coup@(In,113,22,4)
        open(2,file=in,status='old')
        call set_cursor_pos@(1,8)

        WRITE(5,'(/A)')'  SINGLE COLUMN DATA INPUT [Y]    -1-'
        WRITE(5,'(A,$)')'  TWO COLUMN DATA [X,Y]           -2-   => '
        READ(5,*)ISTLP
        IF (ISTLP.EQ.1)THEN
          WRITE(5,'(/A,$)')'  TO ADD 1.COLUMN WITH EQUIDISTANT X-VALUES
     PINTO OUTPUT:  N-NO, Y-YES => '
          READ(5,'(A)')KEY
          IF(KEY.EQ.'y'.OR.KEY.EQ.'Y')THEN
            ISTLP3=2
          ELSE
            ISTLP3=1
          ENDIF
        ELSE
          ISTLP3=2
        ENDIF
        WRITE(5,'(/a)')'  RECTIFICATION (HAND-FITTING)      -1- '
        WRITE(5,*)' EQUIWALENT WIDTH                  -2- '
        WRITE(5,*)' TREND SUBSTRACTION                -3- '
        WRITE(5,'(A,$)')'  RADIAL VELOCITY                   -4-    => '
        READ(5,*)IVOL2
        IF (IVOL2.EQ.1)THEN
          IVOL3=0
        ENDIF
        IF (IVOL2.EQ.2)THEN
          IVOL3=0
        ENDIF
        IF (IVOL2.EQ.3)THEN
          IVOL3=1
          IVOL2=1
        ENDIF
        IF (IVOL2.EQ.4)THEN
          IVOL3=1
          IVOL2=2
        ENDIF
C
C                       IVOL2   IVOL3
C       RECTIF            1       0
C       TREND             1       1
C       EWQ               2       0
C       RV                2       1
C
        IF (IVOL2.EQ.2)THEN
C         write(5,'(/a,$)')'  OUTPUT FILE    => '
C         read(5,'(a12)')out
C         open(3,file=OUT,status='UNKNOWN')
          write(5,'(/a,$)')'  OUTPUT FILE    => '
          read(5,'(a12)')out
106       INQUIRE(FILE=out,EXIST=ex)
          if(.NOT.ex) then
            open(3,file=out,status='new')
          else
            write(5,11)out
            write(5,12)
            read(5,'(a)') key1
            if(key1.eq.'n'.OR.key1.eq.'N') then
              write(5,'(a,$)')' Output file =>'
              read(5,'(a)')out
              goto 106
            else
              open(3,file=out,status='old')
            end if
          end if
        ENDIF
15      IF (IVOL2.EQ.2)THEN
          IF (ISTLP.EQ.2)THEN
            write(5,'(/a,$)')'  REGION FROM-TO     [Xmin,Xmax]  => '
          ENDIF
          IF (ISTLP.EQ.1)THEN
            write(5,'(/a,$)')'  REGION FROM-TO     [Nmin-1,Nmax+1]  => '
          ENDIF
          read(5,*)xxd,xxu
        ELSE
          XXD=-1.D6
          XXU=1.D6
        ENDIF
        N1=0
        N=0
        IF(ISTLP.EQ.1)THEN
          write(5,'(/a,$)')'  START-END  LAMBDA [Lmin,Lmax]  => '
          read(5,*)FREQ1,FREQ2
20        N1=N1+1
          READ(2,*,END=21)YY
          XX=DBLE(N1)
          IF(xx.gt.xxd)then
            IF(XX.LT.XXU)THEN
              N=N+1
              x(n)=xx
              y(n)=yy
            ELSE
              GOTO 21
            ENDIF
          ENDIF
          goto 20
21        DO 22 IFR=1,N
            X(IFR)=(FREQ2-FREQ1)/DBLE(N-1)*DBLE(IFR-1)+FREQ1
22        CONTINUE
          GOTO 30
        ENDIF
        IF(ISTLP.EQ.2)THEN
23        read(2,*,end=30)xx,yy
          if(xx.gt.xxd)then
            IF(XX.LT.XXU)THEN
              N=N+1
              x(n)=xx
              y(n)=yy
            ELSE
              GOTO 30
            ENDIF
          endif
          goto 23
        ENDIF
30      WRITE(5,'(/A,I5)')'  NUMBER OF DATA',N
        WRITE(5,'(10X,A)')'REL. ERROR OF GRAPHICS =3/1000'
        X1=X(1)
        X2=X(1)
        Y1=Y(1)
        Y2=Y(1)
        DO 40 I=1,N
          IF(X(I).LT.X1)X1=X(I)
          IF(X(I).GT.X2)X2=X(I)
          IF(Y(I).LT.Y1)Y1=Y(I)
          IF(Y(I).GT.Y2)Y2=Y(I)
40      CONTINUE
        WRITE(5,'(/A,4E11.3)')'  XMIN XMAX YMIN YMAX',X1,X2,Y1,Y2
        REWIND(2)
C        PAUSE
C       INPUT:N,X1,X2,Y1,Y2
C       OUTPUT:XP,YP,NINTRP SURADNICE FITU A POCET ZADANYCH BODOV FITU
        CALL GRAF(N,XP,YP,NINTRP,X1,X2,Y1,Y2,IVOL2,IVOL3)
C
        DO 107 I=1,NINTRP
          WRITE(4,*)XP(I),YP(I)
107     CONTINUE
        WRITE(4,*)
        IF (IVOL2.EQ.1)THEN
          write(5,'(a,$)')'  OUTPUT FILE    => '
          read(5,'(a12)')out
109       INQUIRE(FILE=out,EXIST=ex)
          if(.NOT.ex) then
            open(3,file=out,status='new')
          else
            write(5,11)out
            write(5,12)
            read(5,'(a)') key1
            if(key1.eq.'n'.OR.key1.eq.'N') then
              write(5,'(a,$)')' Output file =>'
              read(5,'(a)')out
              goto 109
            else
              open(3,file=out,status='old')
            end if
          end if
11        format(5x,'    FILE:    ','[1m',a,'[0m',' exists! ')
12        format(1x,
     + '[1mN[0mew name OR [1mO[0mverwrite? [1m[N/O] [0m (O) => ',
     + $)
          DO 50 I=1,N
            CALL INTEP(X(I),Y0,XP,YP,NINTRP,IER)
            IF (IVOL3.EQ.0)THEN
              Y(I)=Y(I)/Y0
            ELSE
              Y(I)=Y(I)-Y0
            ENDIF
            IF(ISTLP3.EQ.1)THEN
              WRITE(3,*)Y(I)
            ELSE
              WRITE(3,*)X(I),Y(I)
            ENDIF
50        CONTINUE
        ENDIF
C
        WRITE(5,'(/A)')'  YOUR FITTING-POINTS WILL APPEAR IN FILE "FIT"'
        WRITE(5,'(/A)')'  End               -E-'
        IF (IVOL2.EQ.2) WRITE(5,'(A)')'  new Region        -R-'
        WRITE(5,'(A,$)')'  new Input file    -I-  => '
        READ(5,'(A1)')KEY
        IF(KEY.EQ.'I'.OR.KEY.EQ.'i')THEN
          CLOSE(2)
          CLOSE(3)
          GOTO 10
        ENDIF
        IF(KEY.EQ.'R'.OR.KEY.EQ.'r')THEN
          GOTO 15
        ENDIF
        CLOSE(2)
        CLOSE(3)
        CLOSE(4)
        STOP
        END
C-----------------------------------------------------------------------
        SUBROUTINE GRAF(N,XP,YP,NINTRP,XM1,XM2,YM1,YM2,IVOL2,IVOL3)
        PARAMETER (NBOD=20000,NBGR=200)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION IX(NBOD),IY(NBOD),XP(NBGR),YP(NBGR)
        character*12 NAZOV
        LOGICAL EX
        COMMON /A/X(NBOD),Y(NBOD)
        EX=.TRUE.
        IF (IVOL2.EQ.2)THEN
          WRITE(3,'(A,A)')' LINE, EQW, AREA, DELTA LAMBDA  ARE IN YOUR',
     P    ' X-AXIS UNITS'
          WRITE(3,'(A,A)')' LINE                 EQW         AREA',
     P    '    DELTA LAMBDA  VELOCITY [KM/S]'
        ENDIF
        FREQO=1.D5
C       FREQO=0.0D0
        x1=xm1
        x2=xm2
        y1=ym1
        y2=ym2
        XMIN=X1
        XMAX=X2
c       IDLX=bottom left-hand corner x-coordinate, IDLY=y-coordinate
c       IPHX=top right-hand corner x-coordinate, IPHY=y-coordinate
        IDLX=0
        IDLY=430
        IPHX=639
        IPHY=60
        SCAL1=dble(iphX-idlX)/(X2-X1)
        SCAL2=dble(idlY-iphY)/(Y2-Y1)
        SCALX=SCAL1
        NINTRP=0
        KURX=200
        KURY=150
        KROK=49
        CALL VGA@

11      ng=0
        DO 20 J=1,N
        IXs=DINT((X(J)-X1)*SCAL1)+IDLX
        if(ixs.ge.idlX.and.ixs.le.iphX)then
          ng=ng+1
          IX(ng)=ixs
          IY(ng)=-DINT((Y(J)-Y1)*SCAL2)+IDLY
        endif
20      CONTINUE
C
C       PREPISANIE OBRAZOVKY PO KAZDOM ZADANOM BODE
C
150     CALL CLEAR_SCREEN@

        CALL SET_CURSOR_POS@(1,1)
        write(5,'(a,a)')'Home,End,<-->,PgDn,PgUp = to move cursor',
     +  '   +,- = zoom   1-9 = to change step'
        CALL SET_CURSOR_POS@(1,2)
        IF(IVOL2.EQ.1)THEN
        write(5,'(a,a)')'Enter = to set a point',
     +  '   Del = to delete a point   Esc = to finish'
        ELSE
        write(5,'(a,a)')'Enter = to set a point',
     +  '   Esc = to finish'
        ENDIF
C               OS X
        CALL DRAW_LINE@(IDLX,IDLY,IPHX,IDLY,15)
        CALL DRAW_LINE@(IDLX,IPHY,IPHX,IPHY,15)
C               OS Y
        CALL DRAW_LINE@(IDLX,IDLY,IDLX,IPHY,15)
        CALL DRAW_LINE@(IPHX,IDLY,IPHX,IPHY,15)
C               SPEKTRUM
        CALL POLYLINE@(IX,IY,NG,3)

C               CONTINUUM
        IF(IVOL2.EQ.1)THEN
        IF(NINTRP.GE.2)THEN
          DO 160 IKJ=idlX,iphX,2
            XZ=DBLE(IKJ-IDLX)/SCAL1+X1
            CALL INTEP(XZ,YZ,XP,YP,NINTRP,IER)
            IX1=DINT((XZ-X1)*SCAL1)+IDLX
            IY1=-DINT((YZ-Y1)*SCAL2)+IDLY
            IF(IER.EQ.1)CALL SET_PIXEL@(IX1,IY1,13)
160       CONTINUE
        ENDIF
        ENDIF

C               AKTUALNY KURZOR
90      CALL CROSS(KURX,KURY,3,14)
        XPK=DBLE(KURX-IDLX)/SCAL1+X1
        YPK=DBLE(IDLY-KURY)/SCAL2+Y1
        CALL SET_CURSOR_POS@(3,28)
        write(5,'(a,f10.3,3X,A,E11.4)')'lambda =',xpk,'intensity =',ypk
        IF (IVOL2.EQ.2)THEN
          IF(EX)THEN
            CALL SET_CURSOR_POS@(48,28)
            IF (IVOL3.EQ.0)THEN
              write(5,'(a,F12.6)')'E.W. = ',SUM
            ELSE
              write(5,'(a,F12.6)')'R.V. = ',VEL
            ENDIF
          ENDIF
        ENDIF

C               ZADANE BODY PRE KONTINUUM
        DO 140 IKJ=1,NINTRP
          JX1=DINT((XP(IKJ)-X1)*SCAL1)+IDLX
          JY1=-DINT((YP(IKJ)-Y1)*SCAL2)+IDLY
          if(jx1.ge.idlX.and.jx1.le.iphX)CALL CROSS(JX1,JY1,3,12)
140     CONTINUE

C               TAZISKO CIARY V PRIPADE MERANIA RADIALNYCH RYCHLOSTI
        IF(IVOL2.EQ.2.AND.IVOL3.EQ.1)THEN
C         WRITE(5,*)FREQO,X1,SCAL1,IDLX
C         BACHA JFRM ASI NEMOZE BYT ZAPORNE ANI VACSIE AKO POCET PIXLOV
C         MONITORU!!!
          FRMX1=(FREQO-X1)*SCAL1
          IF(FRMX1.GT.(0-IDLX).AND.FRMX1.LT.(IPHX-IDLX))THEN
            IRV=DINT(FRMX1)+IDLX
            CALL DRAW_LINE@(IRV,IDLY,IRV,IPHY,18)
          ENDIF
        ENDIF

C               SPOJNICA 2 POSLEDNYCH BODOV DEFINUJUCICH CIARU
C         BACHA INTEGER NEMOZE BYT ZAPORNE ANI VACSIE AKO POCET PIXLOV
C         MONITORU!!!
        IF(IVOL2.EQ.2)THEN
        IF(EX)THEN
          AJCX1=(XP(NINTRP-1)-X1)*SCAL1
          AJCY1=(YP(NINTRP-1)-Y1)*SCAL2
          AJCX2=(XP(NINTRP)-X1)*SCAL1
          AJCY2=(YP(NINTRP)-Y1)*SCAL2
          IF(AJCX1.ge.(0-IDLX).and.AJCX1.le.(IPHX-IDLX).AND.
     P      AJCX2.GE.(0-IDLX).AND.AJCX2.LE.(IPHX-IDLX))THEN
            JCX1=DINT(AJCX1)+IDLX
            JCY1=-DINT(AJCY1)+IDLY
            JCX2=DINT(AJCX2)+IDLX
            JCY2=-DINT(AJCY2)+IDLY
            CALL DRAW_LINE@(JCX1,JCY1,JCX2,JCY2,17)
          ENDIF
        ENDIF
        ENDIF

100     CALL GET_KEY@(IKEY)
C       RETURN,ESC,->,<-,HORE,DOLE
C       ---------------NACITANIE BODU (ENTER)-------------
        IF(IKEY.EQ.13.AND.NINTRP.LT.NBGR)THEN
          NINTRP=NINTRP+1
          XP(NINTRP)=DBLE(KURX-IDLX)/SCAL1+X1
          YP(NINTRP)=DBLE(IDLY-KURY)/SCAL2+Y1
          IF (IVOL2.EQ.2)THEN
            CALL SET_CURSOR_POS@(48,28)
            EX=.NOT.EX
            IF(EX)THEN
              CALL WIDTH (XP(NINTRP-1),YP(NINTRP-1)
     P        ,XP(NINTRP),YP(NINTRP),N,SUM,FREQO)
              AREA=SUM*(YP(NINTRP-1)+YP(NINTRP))*0.5
              DFREQ=FREQO-FREQL
              VEL=DFREQ/FREQO*2.998D5
              IF (IVOL3.EQ.0)THEN
                WRITE(3,'(A12,3X,F12.6,3X,E11.4,3X,F10.4,3X,F10.1)')
     P          NAZOV,SUM,AREA,DFREQ,VEL
              ELSE
                WRITE(3,'(F12.4,3X,F12.6,3X,E11.4,3X,F10.4,3X,F10.1)')
     P          FREQL,SUM,AREA,DFREQ,VEL
              ENDIF
            ELSE
              IF (IVOL3.EQ.0)THEN
                WRITE(5,'(A,$)')'LINE NAME[12 CHAR]=>'
                READ(5,'(A12)')NAZOV
                FREQL=0.
              ELSE
                WRITE(5,'(A,$)')'LAMBDA LABORATORY =>'
                READ(5,*)FREQL
C                NAZOV=' RADIAL VEL '
              ENDIF
            ENDIF
          ENDIF
          IF (IVOL2.EQ.1)THEN
C               ZNOVUZORAD BODY
         do 210 i=1,nintrp-1
          do 210 j=i+1,nintrp
            if(xp(j).lt.xp(i))then
              a=xp(j)
              xp(j)=xp(i)
              xp(i)=a
              a=yp(j)
              yp(j)=yp(i)
              yp(i)=a
            endif
210       continue
          ENDIF
          KROK=13
          GOTO 150
        ENDIF
C       ---------------ODMAZANIE BODU (Del)-------------
        IF (IVOL2.EQ.1)THEN
        IF(IKEY.EQ.339.AND.NINTRP.GT.1)THEN
          XKUR=DBLE(KURX-IDLX)/SCAL1+X1
          VZDA=1.D32
          DO 350 I=1,NINTRP
            VZ=DABS(XKUR-XP(I))
            IF(VZ.LT.VZDA)THEN
              VZDA=VZ
              IODMAZ=I
            ENDIF
350       CONTINUE
          DO 351 I=IODMAZ,NINTRP-1
            YP(i)=YP(i+1)
            XP(i)=XP(i+1)
351       CONTINUE
          NINTRP=NINTRP-1
          GOTO 150
        ENDIF
        ENDIF
C       ---------------POSUN NA DALSIU STRANU( PgUp )-------------
1       IF(IKEY.EQ.329)THEN
          XPZ=DBLE(IPHX-IDLX)/SCAL1
          IF(X1+XPZ.LE.XMAX-XPZ)X1=X1+XPZ
          GOTO 11
        ENDIF
C       ---------------POSUN NA PREDCHADZAJUCU STRANU( PgDn )-------------
2       IF(IKEY.EQ.337)THEN
          XPZ=DBLE(IPHX-IDLX)/SCAL1
          X1=X1-XPZ
          IF(X1.LT.XMIN)X1=XMIN
          GOTO 11
        ENDIF
C       ---------------POSUN NA ZACIATOK  ( Home )-------------
        IF(IKEY.EQ.327)THEN
          X1=XMIN
          GOTO 11
        ENDIF
C       ---------------POSUN NA KONIEC  ( End )-------------
        IF(IKEY.EQ.335)THEN
          XPZ=DBLE(IPHX-IDLX)/SCAL1
          X1=XMAX-XPZ
          GOTO 11
        ENDIF
C       ---------------ZOOM ZVACSENIE  ( + )-------------
        IF(IKEY.EQ.43)THEN
          IF((X(N)-X(1))*2.D0*SCAL1.GT.32000.D0)THEN
            SCAL1=SCAL1
            GOTO 11
          ELSE
            SCAL1=SCAL1*2.D0
            GOTO 11
          ENDIF
        ENDIF
C       ---------------ZOOM ZMENSENIE  ( - )-------------
        IF(IKEY.EQ.45)THEN
          SCAL1=SCAL1/2.D0
          IF(SCAL1.LT.SCALX)THEN
           SCAL1=SCALX
           X1=XMIN
          ENDIF
          GOTO 11
        ENDIF
C       ---------------POSUN VPRAVO ( -> )-------------
        IF(IKEY.EQ.333)THEN
          CALL CROSS(KURX,KURY,3,0)
          KURX=KURX+KROK
          IF(KURX.GT.IPHX)THEN
           KURX=IDLX
           IKEY=329
           GOTO 1
          ENDIF
          GOTO 90
        ENDIF
C       ---------------POSUN VLAVO ( <- )-------------
        IF(IKEY.EQ.331)THEN
          CALL CROSS(KURX,KURY,3,0)
          KURX=KURX-KROK
          IF(KURX.lT.IDLX)THEN
           KURX=IPHX
           IKEY=337
           GOTO 2
          ENDIF
          GOTO 90
        ENDIF
C       ---------------POSUN HORE (  )-------------
        IF(IKEY.EQ.328)THEN
          CALL CROSS(KURX,KURY,3,0)
          KURY=KURY-KROK
          IF(KURY.LT.IPHY)KURY=IDLY
          GOTO 90
        ENDIF
C       ---------------POSUN DOLE (  )-------------
        IF(IKEY.EQ.336)THEN
          CALL CROSS(KURX,KURY,3,0)
          KURY=KURY+KROK
          IF(KURY.GT.IDLY)KURY=IPHY
          GOTO 90
        ENDIF
C       ---------------ZMENA KROKU ( 1-9 )-------------
        IF(IKEY.GE.49.AND.IKEY.LE.57)THEN
          KROK=(IKEY-48)*6-5
          GOTO 100
        ENDIF
C       ---------------KONIEC ( Esc )----------
        IF(IKEY.NE.27)GOTO 100
C-----------------------------------------------------------------------
        CALL TEXT_MODE@
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE CROSS(IX,IY,ILENG,ICOLL)
        CALL DRAW_LINE@(IX-ILENG,IY,IX+ILENG,IY,ICOLL)
        CALL DRAW_LINE@(IX,IY-ILENG,IX,IY+ILENG,ICOLL)
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE INTEP(XP,P,X,F,N,IER)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       SUBROUTINE   I N T E P  : SPLINE INTERPOLATION SCHEME BASED   C
C                                 ON HERMITE POLYNOMIALS              C
C        INPUT : XP  - THE CHOSEN ARGUMENT VALUE                      C
C                X   - THE VECTOR OF INDEPENDENT VALUES               C
C                F   - THE VECTOR OF FUNCTION OR DEPENDENT VALUES     C
C                N   - THE NUMBER OF POINTS IN THE (X,F) VECTORS      C
C       OUTPUT : P   - THE RESULTANT INTERPOLATED VALUE               C
C                IER - THE RESULTANT ERROR PARAMETER                  C
C       REFERENCE : PUBL.DOM.AST.OBS.VICTORIA,VOL.XVI,NO.6,1982       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        REAL LP1,LP2,L1,L2
        DIMENSION F(N),X(N)
        IER=1
        IO=1
        IUP=0
        IF(X(2).LT.X(1)) IUP=1
        N1=N-1
        IF((XP.GE.X(N).AND.IUP.EQ.0).OR.(XP.LE.X(N).AND.IUP.EQ.1)) THEN
         P=F(N)
         IER=2
         RETURN
        ELSE IF((XP.LE.X(1).AND.IUP.EQ.0).OR.
     *   (XP.GE.X(1).AND.IUP.EQ.1)) THEN
         P=F(1)
         IER=2
         RETURN
        END IF
        DO 1 I=IO,N
         IF(XP.LT.X(I).AND.IUP.EQ.0) GOTO 2
         IF(XP.GT.X(I).AND.IUP.EQ.1) GOTO 2
1       CONTINUE
        P=F(N)
        IER=2
        RETURN
2       I=I-1
        IF(I.EQ.IO-1) GOTO 4
        IO=I+1
        LP1=1./(X(I)-X(I+1))
        LP2=1./(X(I+1)-X(I))
        IF(I.EQ.1) FP1=(F(2)-F(1))/(X(2)-X(1))
        IF(I.EQ.1) GOTO 3
        FP1=(F(I+1)-F(I-1))/(X(I+1)-X(I-1))
3       IF(I.GE.N1) FP2=(F(N)-F(N-1))/(X(N)-X(N-1))
        IF(I.GE.N1) GOTO 4
        FP2=(F(I+2)-F(I))/(X(I+2)-X(I))
4       XPI1=XP-X(I+1)
        XPI=XP-X(I)
        L1=XPI1*LP1
        L2=XPI*LP2
        P=F(I)*(1.-2.*LP1*XPI)*L1*L1+F(I+1)*(1.-2.*LP2*XPI1)
     *  *L2*L2+FP2*XPI1*L2*L2+FP1*XPI*L1*L1
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE WIDTH(XMIN,YMIN,XMAX,YMAX,N,SUM,FREQO)
        PARAMETER (NBOD=20000)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /A/X(NBOD),Y(NBOD)
        SUM=0.D0
        J=0
        DERM=(YMAX-YMIN)/(XMAX-XMIN)
        DO 10 I=1,N
        IF(X(I).GT.XMIN)THEN
          IF(X(I).LT.XMAX)THEN
            J=J+1
C           LOWER EDGE
            IF(J.EQ.1)THEN
              DER=(Y(I)-Y(I-1))/(X(I)-X(I-1))
              DX=X(I)-XMIN
              DSUM=DX*0.5D0*(Y(I)+Y(I-1)+DER*(XMIN-X(I-1)))
              SUM=DSUM
              XMEAN=(X(I)+XMIN)*0.5D0
              RDM=XMEAN*(DSUM-DX*(YMIN+DERM*(XMEAN-XMIN)))
              FREQO=RDM
              GOTO 10
            ENDIF
            DX=X(I)-X(I-1)
            DSUM=DX*(Y(I)+Y(I-1))*0.5D0
            SUM=DSUM+SUM
            XMEAN=(X(I)+X(I-1))*0.5D0
            RDM=XMEAN*(DSUM-DX*(YMIN+DERM*(XMEAN-XMIN)))
            FREQO=RDM+FREQO
          ELSE
            GOTO 20
          ENDIF
        ENDIF
10      CONTINUE
C       UPPER EDGE
20      DER=(Y(I)-Y(I-1))/(X(I)-X(I-1))
        DX=XMAX-X(I-1)
        DSUM=DX*0.5D0*(2.D0*Y(I-1)+DER*DX)
        SUM=DSUM+SUM
        XMEAN=(X(I-1)+XMAX)*0.5D0
        RDM=XMEAN*(DSUM-DX*(YMIN+DERM*(XMEAN-XMIN)))
        FREQO=RDM+FREQO
C         SUBSTRACT THE CONTINUUM AREA AND GET THE LINE AREA
        SUM=(XMAX-XMIN)*0.5D0*(YMAX+YMIN)-SUM
        FREQO=FREQO/SUM
        IF(FREQO.LT.0.)FREQO=-FREQO
C         GET EQUIVALENT WIDTH
        SUM=SUM/(YMAX+YMIN)*2.D0
C         SUM IS EQUIVALENT WIDTH OF THE SPECTRAL LINE
C         SUM>0 FOR ABSORPTION, SUM<0 FOR EMISSION
C         FREQO IS WAVELENGHT OF THE MASS CENTER OF THE LINE
        RETURN
        END
