C       OPTIONS(SAVE,DREAL,CHECK)
       OPTIONS(SAVE,DREAL,OPTIMISE)
C       THESE ARE OPTIONS FOR THE FTN77 COMPILER UNIVERSITY OF SALFORD
C       REMOVE THEM USING C IN THE FIRST COLUMN (COMMENT) IF USING
C       OTHER COMPILER.
        PROGRAM KONVOL7
C               ~~~~~~~
C       VERSION: JULY 3, 1995
C       PROGRAMM FOR CALCULATING OF THE CONVOLUTION WITH CHOSEN FUNCTION
C       BASED ON NUMERICAL INTEGRATION.
C       FOR YOUR OWN FILTER FUNCTION JUST WRITE YOUR OWN
C                        'FUNCTION FILTER(X)'
C       RENAME THE ORIGINAL ONE AND COMPILE OR USE THE OPTION
C                       'OWN TABULATED FILTER'.
C       YOU ARE SUPPOSED TO CHOSE AREA OF YOUR OWN TABULATED FILTER=1
C       NOT TO CHANGE AN INTEGRAL OF THE ORIGINAL FUNCTION
C       AND TO OPERATE GRAPHICS.
C       ALL BUILD IN FILTERS ARE AUTOMATICALLY NOTMALIZED.
c
c       authors: J. Budaj
c       contact: budaj@ta3.sk
c       published: http://www.ta3.sk/~budaj/software
c       special requirements: DBOS environement,
c       non original routins:  intep, hunt, hunt1 (Press.et.al. 1992,
c                                           Numerical recipes)
c                                LIMBDARK (WRITTEN BY BARONE ET AL.)
c                                fint (WRITTEN BY BARONE ET AL.)
c       input: keyboard + files
c       output:  file
C
C       REMARKS: ROTATIONAL HALFWIDTH CHANGES AT EACH POINT
C       OPTION EXPLANATION:
C         SLOW,LITTLE DATA IN FILTER WIDTH <2000     -1- :
C               PRECISE-HERMIT POLYNOMES INTERPOLATION, DATA WITH GAPES
C         FAST,UNIVERSAL      <20000                 -2- :
C               LINEAR INTERPOLATION, DATA WITH GAPES
C         FASTER,ENOUGH DATA IN FILTER WIDTH <20000  -3- :
C               LINEAR INTERPOLATION, NOT OPTIONAL STEP=STEPS IN DATA
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        CHARACTER*12 IN,INN,INP
        CHARACTER*30 CURDIR@,IPA
        COMMON /A/ X1(20000),F1(20000)
        COMMON /IVO/ IVOL
        COMMON /CONST/SPACE,NSPACE,CUTOFF,NDAT1
        COMMON /ALF1/ A,AA,AAA
        COMMON /ALF2/ DLL,UKOEF
        COMMON /BET/ FLTRX(1000),FLTRY(1000),NDAT2
        COMMON /B/SCAL1,SCAL2,IDL1,IDL2,XMIN,YMIN,XMAX,YMAX
C       REMOVE C IN THE FIRST COLUMN (COMMENT) IN THE NEXT LINE
C       IF USING OTHER COMPILER
C       OPEN(5,FILE='CON')
C       PUT 'C' INTO THE FIRST COLUM OF THE THE NEXT 13 LINES
C       IF USING OTHER COMPILER
        CALL CLEAR_SCREEN@
        CALL HIDE_CURSOR@
        CALL COUP@(' W E L C O M E   T O ',116,39,20)
        CALL COUP@(' K O N V O L 7 ',244,41,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)
        CALL COUP@(' INPUT FILE NAME : => ',30,0,4)
16      CALL SELECT_FILE@('*.*',IN,*16)
        CALL COUP@(IN,113,22,4)
        CALL SET_CURSOR_POS@(1,8)
        WRITE(5,*)
        WRITE(5,*)'-------------------------------------------'
        WRITE(5,*)'                   PROGRAMM FOR CONVOLUTION'
        WRITE(5,*)'-------------------------------------------'
        WRITE(5,*)
C       REMOVE 'C' IN THE NEXT TWO LINES IF USING OTHER COMPILER
C       WRITE(5,'(A17,$)')'  INPUT FILE   => '
C       READ(5,'(A12)')IN
        WRITE(5,'(A17,$)')'  OUTPUT FILE  => '
        READ(5,'(A12)')INN
5       OPEN(2,FILE=IN,STATUS='OLD')
        OPEN(3,FILE=INN)
C       PUT 'C' IN THE NEXT LINE IF USING OTHER COMPILER
        WRITE(5,*)' TIME',TIME2-TIME1,'SEC'
        WRITE(5,*)
        WRITE(5,*)' GAUSS. FILTER     -1-'
        WRITE(5,*)' ROTATION          -2-'
        WRITE(5,*)' RECTANGLE         -3-'
        WRITE(5,*)' TRIANGLE          -4-'
        WRITE(5,*)' OWN FILTER        -5-'
        WRITE(5,'(A,$)')'  FINISH            -0-      => '
        READ(5,*)IVOL
        IF (IVOL.EQ.0)GOTO 40
        WRITE(5,*)
        WRITE(5,*)
        WRITE(5,*)' NOTE:     WE RECOMMEND TO REACH REL. ERROR 1/1000 '
        IF (IVOL.EQ.1)THEN
        WRITE(5,*)' GAUSS.FILTER: FILTER WIDTH/3=FWHM>1.5*STEP'
        WRITE(5,*)'               FILTER WIDTH/5=EXP(-1)WIDTH>0.9*STEP'
        WRITE(5,*)' FWHM MEANS FULL WIDTH AT HALF MAXIMA'
        ENDIF
        IF (IVOL.EQ.2)THEN
        WRITE(5,*)' ROTATION: STEP<HALFWIDTH/30 (/10?)'
        WRITE(5,*)'           STEP<HALFWIDTH/3  TO REACH REL.ERR. 1/100'
        ENDIF
        IF (IVOL.EQ.3)THEN
        WRITE(5,*)' RECTANGLE: STEP < WIDTH/10'
        ENDIF
        IF (IVOL.EQ.4)THEN
        WRITE(5,*)' TRIANGLE : STEP < BASELINE/30'
        ENDIF
        WRITE(5,*)' BUT ALWAYS: STEP<CHARACTER. DISTANCE'
        WRITE(5,*)' OF THE CHANGES IN YOUR FUNCTION OR FILTER'
        WRITE(5,*)'       [YU] -YOUR UNIT OF THE X AXIS'
        WRITE(5,*)
C        WRITE(5,'(/A)')' SLOW,LITTLE DATA IN FILTER WIDTH <2000     -1-'
C        WRITE(5,'(A)')' FAST,UNIVERSAL      <20000                 -2-'
C        WRITE(5,'(A,$)')' FASTER,ENOUGH DATA IN FILTER WIDTH <20000  -3-
C     P => '
C        READ(5,*)IVOL1
        IVOL1=2
        WRITE(5,*)
        I=1
10      READ(2,*,END=20,ERR=11)X1(I),F1(I)
        I=I+1
        GOTO 10
11      WRITE(5,*)'    WARNING:  INPUT ERROR'
20      NDAT1=I-1
        WRITE(5,'(A,I6)')' No. OF DATA LOADED=',NDAT1
        IF (IVOL.EQ.1)THEN
          WRITE(5,'(/A)')'  FWHM                  -0-'
          WRITE(5,'(A,$)')'  EXP(-1) WIDTH         -1-      => '
          READ(5,*)IVOL2
          IF (IVOL2.EQ.0)THEN
            WRITE(5,'(A,$)')'  FWHM OF THE FILTER [YU]           => '
            READ(5,*)A
C               WE WILL USE F(X)=EXP(-A**2*X**2)
            A=2.D0*DSQRT(-DLOG(0.5D0))/A
            WRITE(5,'(15X,A,E10.3)')' EXP(-1)WIDTH OF THE FILTER',1.D0/A
          ELSE
            WRITE(5,'(A,$)')'  EXP(-1)WIDTH OF THE FILTER [YU]      => '
            READ(5,*)A
            A=1.D0/A
            WRITE(5,'(15X,A,E10.3)')'  FWHM OF THE FILTER',
     P      2.D0*DSQRT(-DLOG(0.5D0))/A
          ENDIF
          WRITE(5,'(A,$)')'  WIDTH OF FILTER [YU]  (AUTOMATIC -0-)  => '
          READ(5,*)CUTOFF
          CUTOFF=CUTOFF*0.5D0
C               IT MEANS WIDTH/3.5=FWHM
          IF(DABS(CUTOFF-0.D0).LT.1.D-30)CUTOFF=2.9D0/A
        ENDIF
        IF (IVOL.EQ.2)THEN
        WRITE(5,'(/A,$)')'  V*SIN I [KM/S], TEMPER.[K], LAMBDA [A]  => '
          READ(5,*)VSIN,TEP,FREQ
C        WRITE(5,'(A,$)')'  WHAT PART OF 1ANGSTREM IS 1YU            => '
C          READ(5,*)UMERA
          CALL LIMBDARK(TEP/10000.D0,FREQ/10000.D0,UKOEF)
C         NOW FREQ IN ANGSTREM, DLL IN YU
          POM1=VSIN/299792.D0
          DLL=POM1*X1(1)
          CUTOFF=DLL
          WRITE(5,*)
          WRITE(5,'(10X,A,F6.0,A,F6.3)')' COEF. OF THE DARKNESS AT',
     P    FREQ,'[A] =',UKOEF
          WRITE(5,'(10X,A,F6.0,A,F8.3)')' ROTATIONAL HALFWIDTH [YU] AT',
     P    X1(1),'[YU] =',DLL
        ENDIF
        IF (IVOL.EQ.3)THEN
          WRITE(5,'(A,$)')'  WIDTH OF THE RECTANGLE [YU]  => '
          READ(5,*)AA
          CUTOFF=0.5D0*AA
          WRITE(5,*)'CUTOFF',CUTOFF
        ENDIF
        IF (IVOL.EQ.4)THEN
          WRITE(5,'(A,$)')'  BASELINE OF THE TRIANGLE [YU] => '
          READ(5,*)AAA
          CUTOFF=0.5D0*AAA
          WRITE(5,*)'CUTOFF',CUTOFF
        ENDIF
        IF (IVOL.EQ.5)THEN
          WRITE(5,'(A,$)')'  FILE NAME WITH FILTER     X,Y DATA   => '
          READ(5,'(A12)')INP
          OPEN(4,FILE=INP,STATUS='OLD')
          I=1
13        READ(4,*,END=14)FLTRX(I),FLTRY(I)
          I=I+1
          GOTO 13
14        NDAT2=I-1
          CLOSE(4)
          CUTOFF=MAX(DABS(FLTRX(NDAT2)),DABS(FLTRX(1)))
          WRITE(5,*)'CUTOFF',CUTOFF
        ENDIF
        IF (IVOL1.NE.3)THEN
          WRITE(5,'(/A,$)')'  STEP [YU]                  => '
          READ(5,*)SPACE0
C               VYROBI NEPARNY POCET BODOV CIZE: NSPACE IS EVEN
          NSPACE=DINT(2.D0*CUTOFF/SPACE0)+1
          IF (DBLE(NSPACE)/2.D0-DBLE(NSPACE/2).GT.0.1) NSPACE=NSPACE+1
          SPACE=2.D0*CUTOFF/DBLE(NSPACE)
          WRITE(5,'(/15X,A,A,A,E11.3,I5,E11.3)')'FILTER WIDTH, ',
     P    ' No. OF STEPS IN IT, ',
     P    ' MODIF.STEP',2.D0*CUTOFF,NSPACE,SPACE
        ENDIF
C                MAX AND MIN VALUES
        XMIN=X1(1)
        XMAX=X1(NDAT1)
        YMIN=F1(1)
        YMAX=F1(1)
        DO 21 I=1,NDAT1
          IF(F1(I).LT.YMIN)YMIN=F1(I)
          IF(F1(I).GT.YMAX)YMAX=F1(I)
21      CONTINUE
        WRITE(5,'(A,A,$)')'  I CACULATE CONVOL. FOR EACH N-TH POINT ',
     P  '[1=FOR EACH] => '
        READ(5,*)NSTEP
C       PUT 'C' IN THE FIRST COLUMN OF THE NEXT TWO LINES IF OTHER COMPILER
        CALL GRAF(NDAT1)
        CALL DCLOCK(TIME1)
        JLO1=1
        JLO2=1
        DO 30 I=1,NDAT1,NSTEP
C         REMARKS: ROTATIONAL HALFWIDTH CHANGES AT EACH POINT
          IF (IVOL.EQ.2)THEN
            DLL=POM1*X1(I)
            CUTOFF=DLL
            NSPACE=DINT(2.D0*CUTOFF/SPACE0)+1
            SPACE=2.D0*CUTOFF/DBLE(NSPACE)
          ENDIF
          IF(X1(I)-X1(1).LT.CUTOFF.OR.X1(NDAT1)-X1(I).LT.CUTOFF)GOTO 30
          IF (IVOL1.EQ.2)THEN
            AK=AKONV2(X1(I),JLO1)
            GOTO 32
          ENDIF
          IF (IVOL1.EQ.1)THEN
            AK=AKONV1(X1(I))
            GOTO 32
          ENDIF
          AK=AKONV3(X1(I),JLO1,JLO2)
32        CONTINUE
C         PUT 'C' IN THE FIRST COLUMN OF THE NEXT LINE IF OTHER COMPILER
          CALL GRAF1(X1(I),AK)
C         WRITE(5,'(E16.8,E14.6)')X1(I),AK
          WRITE(3,'(E16.8,E14.6)')X1(I),AK
30      CONTINUE
C         PUT 'C' IN THE FIRST COLUMN OF THE NEXT LINE IF OTHER COMPILER
        CALL DCLOCK(TIME2)
        CLOSE (2)
        CLOSE (3)
C       PUT 'C' IN THE FIRST COLUMN OF NEXT 2 LINES IF OTHER COMPILER
        CALL GET_KEY@(IVOL3)
        IF(IVOL3.EQ.13)THEN
          GOTO 5
C       PUT 'C' IN THE FIRST COLUMN OF THE NEXT LINE IF OTHER COMPILER
        ENDIF
40      CONTINUE
C        CLOSE(5)
        STOP
        END
C-----------------------------------------------------------------------
        FUNCTION AKONV1(XX)
C       INTERPOLATION OF FUNCTION BY MEANS OF HERMIT POLYNOMS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       X(I)  - POLE X-OVYCH HODNOT, X(I+1) > X(I)
C       F(I)  - POLE Y-OVYCH HODNOT FILTROVANEJ FUNKCIE V BODOCH X(I)
C       XX    - X-OVA SURADNICA BODU KDE TA ZAUJIMA KONVOLUCIA
C       SPACE - INTEGRACNY KROK
C       NSPACE- POCET KROKOV V SIRKE FILTRA
C       CUTOFF- POLOVICA SIRKY FILTRA
C       NDAT1 - POCET DAT
C       AKONV - HODNOTA KONVOLUCIE (VYFITROVANA HODNOTA FUNKCIE) V XX
C       FILTER- FILTRUJUCA FUNKCIA (JE VYHODNE ABY JEJ INTEGRAL=1
C               => FILTROVANIM SA NEMENI INTEGRAL FILTROVANEJ FUNKCIE)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /CONST/SPACE,NSPACE,CUTOFF,NDAT1
        COMMON /OWNFIL/ JLO3
        JLO3=1
        SUMA=0.D0
        XI1=XX-CUTOFF
        CALL INTEP(XI1,YI1,NDAT1,IER)
        F1=FILTER(XI1-XX)
        DO 100 I=1,NSPACE
          XI2=XI1+SPACE/2.D0
          XI3=XI1+SPACE
          F2=FILTER(XI2-XX)
          F3=FILTER(XI3-XX)
          CALL INTEP(XI2,YI2,NDAT1,IER)
          CALL INTEP(XI3,YI3,NDAT1,IER)
          SUMA=(F1*YI1+4.D0*F2*YI2+F3*YI3)+SUMA
          YI1=YI3
          XI1=XI3
          F1=F3
100     CONTINUE
        AKONV1=SPACE*SUMA/6.D0
        RETURN
        END
C-----------------------------------------------------------------------
        FUNCTION AKONV2(XX,JLO1)
C       LINEAR INTERPOLATION OF FUNCTION
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /CONST/SPACE,NSPACE,CUTOFF,NDAT1
        COMMON /OWNFIL/ JLO3
        JLO3=1
        SUMA=0.D0
        XI1=XX-CUTOFF
        CALL HUNT(NDAT1,XI1,JLO1)
        JLO=JLO1
        CALL INTP(JLO,NDAT1,XI1,YI1)
        F1=FILTER(XI1-XX)
        DO 100 I=1,NSPACE
          XI2=XI1+SPACE/2.D0
          XI3=XI1+SPACE
          CALL HUNT(NDAT1,XI2,JLO)
          CALL INTP(JLO,NDAT1,XI2,YI2)
          CALL HUNT(NDAT1,XI3,JLO)
          CALL INTP(JLO,NDAT1,XI3,YI3)
          F2=FILTER(XI2-XX)
          F3=FILTER(XI3-XX)
          SUMA=(F1*YI1+4.D0*F2*YI2+F3*YI3)+SUMA
          YI1=YI3
          XI1=XI3
          F1=F3
100     CONTINUE
        AKONV2=SPACE*SUMA/6.D0
        RETURN
        END
C-----------------------------------------------------------------------
        FUNCTION AKONV3(XX,JLO1,JLO2)
C       VERY FAST INTEGRATION FOR BROAD FILTERS-WITH ENOUGH DATA POINTS
C       IN THE WIDTH OF FILTER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       X1(I)  - POLE X-OVYCH HODNOT, X1(I+1) > X1(I)
C       F1(I)  - POLE Y-OVYCH HODNOT FILTROVANEJ FUNKCIE V BODOCH X1(I)
C       XX    - X-OVA SURADNICA BODU KDE TA ZAUJIMA KONVOLUCIA
C       NDAT1 - POCET DAT
C       AKONV3- HODNOTA KONVOLUCIE (VYFITROVANA HODNOTA FUNKCIE) V XX
C       FILTER- FILTRUJUCA FUNKCIA (JE VYHODNE ABY JEJ INTEGRAL=1
C               => FILTROVANIM SA NEMENI INTEGRAL FILTROVANEJ FUNKCIE)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C        COMMON /A/ X(20000),F(20000)
        COMMON /A/ X1(20000),F1(20000)
        COMMON /CONST/SPACE,NSPACE,CUTOFF,NDAT1
        COMMON /OWNFIL/ JLO3
        JLO3=1
        XI1=XX-CUTOFF
        XI2=XX+CUTOFF
        CALL HUNT(NDAT1,XI1,JLO1)
        CALL HUNT(NDAT1,XI2,JLO2)
        IF((JLO2-JLO1).LT.2)THEN
          WRITE(5,'(A)')'  NOT ENOUGH DATA IN THE FILTER WIDTH'
          STOP
        ENDIF
        SUMA=0.D0
        SUMA=F1(JLO1+1)*FILTER(X1(JLO1+1)-XX)*0.5D0
     P   *(X1(JLO1+1)-XI1)+SUMA
        SUMA=F1(JLO2)*FILTER(X1(JLO2)-XX)*0.5D0
     P   *(XI2-X1(JLO2))+SUMA
        DO 100 I=JLO1+1,JLO2-1
          F10=F1(I)
          F20=F1(I+1)
          SUMA=(F20*FILTER(X1(I+1)-XX)+F10*FILTER(X1(I)-XX))*0.5D0
     P    *(X1(I+1)-X1(I))+SUMA
100     CONTINUE
        AKONV3=SUMA
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE INTP(JLO,NDAT1,XI,YI)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /A/ X1(20000),F1(20000)
C       LINEAR INTERPOLATION
        YI=(F1(JLO+1)-F1(JLO))*(XI-X1(JLO))/(X1(JLO+1)-X1(JLO))+F1(JLO)
        RETURN
        END
C-----------------------------------------------------------------------
        FUNCTION FILTER(X)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /IVO/ IVOL
        IF (IVOL.EQ.1)THEN
          FILTER=GAUS(X)
          RETURN
        ENDIF
        IF (IVOL.EQ.2)THEN
          FILTER=ROT(X)
          RETURN
        ENDIF
        IF (IVOL.EQ.3)THEN
          FILTER=RECT(X)
          RETURN
        ENDIF
        IF (IVOL.EQ.4)THEN
          FILTER=TRIAN(X)
          RETURN
        ENDIF
        IF (IVOL.EQ.5)THEN
          FILTER=OWN(X)
          RETURN
        ENDIF
        END
C-----------------------------------------------------------------------
        FUNCTION OWN(X)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /BET/ FLTRX(1000),FLTRY(1000),NDAT2
        COMMON /OWNFIL/ JLO3
        IF(X.LT.FLTRX(1).OR.X.GT.FLTRX(NDAT2))THEN
          OWN=0.D0
          RETURN
        ENDIF
        CALL HUNT1(NDAT2,X,JLO3)
        OWN=(X-FLTRX(JLO3))*(FLTRY(JLO3+1)-FLTRY(JLO3))/
     P  (FLTRX(JLO3+1)-FLTRX(JLO3))+FLTRY(JLO3)
        RETURN
        END
C-----------------------------------------------------------------------
        FUNCTION GAUS(X)
C               1/A IS WIDTH AT A/SQRT(3.14)*EXP(-1)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /ALF1/ A,AA,AAA
        AL=A*A*X*X
        IF(AL.GT.16.D0)THEN
          GAUS=0.D0
          RETURN
        ELSE
C               1/DSQRT(3.1415927) =0.5641896
          GAUS=A*0.5641896D0*DEXP(-AL)
          RETURN
        ENDIF
        END
C-----------------------------------------------------------------------
        FUNCTION RECT(X)
C               AA IS THE WIDTH OF THE RECTANGLE
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /ALF1/ A,AA,AAA
        IF(DABS(X).GT.AA)THEN
          RECT=0.D0
          RETURN
        ELSE
C               NORMALIZATION
          RECT=1.0D0/AA
          RETURN
        ENDIF
        END
C-----------------------------------------------------------------------
        FUNCTION TRIAN(X)
C               AAA IS THE BASELINE OF THE TRIANGLE
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /ALF1/ A,AA,AAA
C               NORMALIZATION
        HEIGHT=2.D0/AAA
        ABX=DABS(X)
        IF(ABX.GT.AAA)THEN
          TRIAN=0.D0
          RETURN
        ELSE
          TRIAN=HEIGHT-2.D0*HEIGHT/AAA*ABX
          RETURN
        ENDIF
        END
C-----------------------------------------------------------------------
        FUNCTION ROT(DL)
C               DL   - ARGUMENT ROTACNEJ FUNKCIE V [TJ]
C               DLL  - ROTACNA POLOSIRKA V [TJ]
C               UKOEF- KOEF.STEMNENIA
C       VYSTUP       - ROT. PROFIL V [TJ]
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /ALF2/ DLL,UKOEF
        IF(DL*DL.GT.DLL*DLL)THEN
          ROT=0.D0
          RETURN
        ENDIF
C        ROT=(2.D0*(1.D0-UKOEF)*DSQRT(1.D0-DL*DL/DLL/DLL)+0.5D0*3.1415D0
C     P  *UKOEF*(1.D0-DL*DL/DLL/DLL))/(3.1415D0*DLL*(1.D0-UKOEF/3.D0))
        ROT=(0.6366198D0*(1.D0-UKOEF)*DSQRT(1.D0-DL*DL/DLL/DLL)+0.5D0
     P  *UKOEF*(1.D0-DL*DL/DLL/DLL))/(DLL*(1.D0-UKOEF*0.3333333D0))
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE LIMBDARK(TEMP,WAVLEN,XX)
C THIS SUBROUTINE WAS WRITTEN BY BARONE ET AL.
C
C TEMP   = TEMPERATURE, IN UNITS OF 10000 K
C WAVLEN = WAVELENGTH, IN UNITS OF ANGSTROM * 10000
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION ARG(2),NENT(2),ENT(35),TABLE(294)
     P  ,TABLE1(105),TABLE2(105),TABLE3(84)
C
C       THE TABLE IS FORMATTED AS FOLLOWS:
C
C --------------------------------------
C            T1   T2   T3   T4   T5  ...
C --------------------------------------
C    WL1      1    2    3    4    5
C    WL2      6    7    8    9   10
C    WL3     11   12   13   14   15
C    ...
C    ...
C --------------------------------------
C
C
        DATA TABLE1
     1  /0.06,0.18,0.52,0.52,0.28,1.00,1.00,1.00,1.00,1.00,
     2  1.00,1.00,0.98,0.93,0.87,0.79,0.73,0.68,0.65,0.58,0.51,
     3  0.99,0.99,1.00,0.97,0.90,0.82,0.78,0.74,0.70,0.64,
     4  0.60,0.59,0.57,0.56,0.53,0.49,0.46,0.44,0.43,0.40,0.38,
     5  0.87,0.87,0.95,0.84,0.76,0.68,0.63,0.58,0.53,0.47,
     6  0.42,0.40,0.39,0.40,0.38,0.36,0.35,0.34,0.34,0.33,0.32,
     7  0.97,0.97,0.94,0.87,0.83,0.80,0.78,0.74,0.72,0.69,
     8  0.63,0.59,0.55,0.52,0.49,0.45,0.43,0.41,0.39,0.37,0.34,
     9  1.00,0.99,0.85,0.78,0.75,0.72,0.71,0.68,0.66,0.62,
     1  0.57,0.54,0.51,0.48,0.47,0.42,0.40,0.38,0.36,0.34,0.32/
       DATA TABLE2
     1 /0.97,0.90,0.77,0.71,0.68,0.66,0.65,0.63,0.60,0.56,
     3 0.52,0.49,0.47,0.44,0.41,0.38,0.35,0.34,0.32,0.31,0.29,
     4 0.88,0.83,0.71,0.65,0.62,0.60,0.59,0.57,0.55,0.51,
     5 0.47,0.45,0.43,0.40,0.38,0.34,0.32,0.31,0.30,0.28,0.27,
     6 0.81,0.76,0.65,0.60,0.57,0.55,0.54,0.52,0.50,0.46,
     7 0.43,0.41,0.39,0.37,0.35,0.31,0.30,0.28,0.27,0.26,0.25,
     8 0.69,0.65,0.56,0.52,0.49,0.47,0.46,0.44,0.42,0.38,
     9 0.36,0.34,0.33,0.31,0.29,0.27,0.25,0.24,0.23,0.22,0.22,
     1 0.61,0.58,0.50,0.49,0.43,0.41,0.40,0.38,0.36,0.33,
     2 0.30,0.29,0.28,0.26,0.25,0.23,0.21,0.20,0.20,0.20,0.20/
       DATA TABLE3
     3 /0.52,0.49,0.43,0.40,0.38,0.37,0.37,0.36,0.34,0.31,
     4 0.29,0.28,0.26,0.25,0.24,0.21,0.20,0.19,0.18,0.18,0.18,
     5 0.48,0.46,0.40,0.36,0.35,0.33,0.32,0.31,0.28,0.26,
     6 0.24,0.23,0.22,0.21,0.20,0.18,0.17,0.16,0.15,0.15,0.16,
     7 0.40,0.38,0.33,0.31,0.29,0.28,0.26,0.25,0.23,0.21,
     8 0.20,0.19,0.18,0.17,0.16,0.14,0.13,0.13,0.12,0.12,0.13,
     9 0.33,0.31,0.27,0.25,0.23,0.22,0.20,0.19,0.17,0.16,
     1 0.14,0.13,0.13,0.12,0.11,0.10,0.09,0.09,0.09,0.10,0.10/
        DATA NENT/21,14/
        DATA ENT/4000.,4500.,5000.,5500.,6000.,6500.,7000.,7500.,
     1 8000.,8500.,9000.,9500.,10000.,11000.,12000.,14000.,
     2 16000.,18000.,20000.,25000.,30000.,
     3 2000.,3000.,3600.,4000.,4500.,5000.,5500.,6000.,7000.,
     4 8000.,10000.,12000.,16000.,22000./
        DO 10 I=1,294
          IF(I.LE.105)TABLE(I)=TABLE1(I)
          IF(I.GE.106.AND.I.LE.210)TABLE(I)=TABLE2(I-105)
          IF(I.GE.211)TABLE(I)=TABLE3(I-210)
10      CONTINUE
        NARG=2
        ARG(1)=TEMP*10000.
        ARG(2)=WAVLEN*10000.
        XX=FINT(NARG,ARG,NENT,ENT,TABLE)
        RETURN
        END
C-----------------------------------------------------------------------
       FUNCTION FINT(NARG,ARG,NENT,ENT,TABLE)
C THIS SUBROUTINE WAS WRITTEN BY BARONE ET AL.
C**     FUNCTION PER INTERPOL. LINEARE IN UNA TABELLA MULTIDIMENSIONALE
C
C       CHIAMATA:       VALINT= FINT(NARG,ARG,NENT,ENT,TABLE)
C
C       ARGOMENTI:      NARG= DIMENSIONE DEL VETTORE ARG (<6)
C                       ARG= COORDINATE DEL PUNTO IN CUI INTERPOLARE
C                       NENT= LISTA DEI NUMERI DI ENTRIES PER OGNI ARG.
C                       ENT= ENTRIES DELLA TABELLA
C                       TABLE= VALORI DELLA TABELLA
C
C      REAL*8 FUNCTION FINT(NARG,ARG,NENT,ENT,TABLE)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ARG(2),NENT(2),ENT(35),TABLE(294)
C      DIMENSION ARG(5),NENT(5),ENT(10),TABLE(10)   PRE F77 VER3.34
      DIMENSION D(5),NCOMB(5),IENT(5)
      KD=1
      M=1
      JA=1
      DO 80 I=1,NARG
      NCOMB(I)=1
      JB=JA-1+NENT(I)
      DO 20 J=JA,JB
      IF (ARG(I).LE.ENT(J)) GO TO 40
20    CONTINUE
      J=JB
40    IF (J.NE.JA) GO TO 60
      J=J+1
60    JR=J-1
      D(I)=(ENT(J)-ARG(I))/(ENT(J)-ENT(JR))
      IENT(I)=J-JA
      KD=KD+IENT(I)*M
      M=M*NENT(I)
80    JA=JB+1
      FINT=0.
100   FAC=1.
      IADR=KD
      IFADR=1
      DO 140 I=1,NARG
      IF (NCOMB(I).EQ.0) GO TO 120
      FAC=FAC*(1.-D(I))
      GO TO 140
120   FAC=FAC*D(I)
      IADR=IADR-IFADR
140   IFADR=IFADR*NENT(I)
      FINT=FINT+FAC*TABLE(IADR)
      IL=NARG
160   IF (NCOMB(IL).EQ.0) GO TO 200
      NCOMB(IL)=0
      IF (IL.EQ.NARG) GO TO 100
      IL=IL+1
      DO 180  K=IL,NARG
180   NCOMB(K)=1
      GO TO 100
200   IL=IL-1
      IF(IL.NE.0) GO TO 160
      RETURN
      END
C-----------------------------------------------------------------------
C       SUBROUTINE   I N T E P  : SPLINE INTERPOLATION SCHEME BASED
C                                 ON HERMITE POLYNOMIALS
C        INPUT : XP  - THE CHOSEN ARGUMENT VALUE
C                X   - THE VECTOR OF INDEPENDENT VALUES
C                F   - THE VECTOR OF FUNCTION OR DEPENDENT VALUES
C                N   - THE NUMBER OF POINTS IN THE (X,F) VECTORS
C       OUTPUT : P   - THE RESULTANT INTERPOLATED VALUE
C                IER - THE RESULTANT ERROR PARAMETER
        SUBROUTINE INTEP(XP,P,N,IER)
C       SUBROUTINE INTEP(XP,P,X,F,N,IER)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        REAL LP1,LP2,L1,L2
C       DIMENSION F(NDIM),X(NDIM)
C       DIMENSION F(N),X(N)
        COMMON /A/ X1(20000),F1(20000)
        IER=1
        IO=1
        IUP=0
        IF(X1(2).LT.X1(1)) IUP=1
        N1=N-1
        IF((XP.GE.X1(N).AND.IUP.EQ.0).OR.(XP.LE.X1(N).AND.IUP.EQ.1))THEN
         P=F1(N)
         IER=2
         RETURN
        ELSE IF((XP.LE.X1(1).AND.IUP.EQ.0).OR.
     *   (XP.GE.X1(1).AND.IUP.EQ.1)) THEN
         P=F1(1)
         IER=2
         RETURN
        END IF
        DO 1 I=IO,N
         IF(XP.LT.X1(I).AND.IUP.EQ.0) GOTO 2
         IF(XP.GT.X1(I).AND.IUP.EQ.1) GOTO 2
1       CONTINUE
        P=F1(N)
        IER=2
        RETURN
2       I=I-1
        IF(I.EQ.IO-1) GOTO 4
        IO=I+1
        LP1=1./(X1(I)-X1(I+1))
        LP2=1./(X1(I+1)-X1(I))
        IF(I.EQ.1) FP1=(F1(2)-F1(1))/(X1(2)-X1(1))
        IF(I.EQ.1) GOTO 3
        FP1=(F1(I+1)-F1(I-1))/(X1(I+1)-X1(I-1))
3       IF(I.GE.N1) FP2=(F1(N)-F1(N-1))/(X1(N)-X1(N-1))
        IF(I.GE.N1) GOTO 4
        FP2=(F1(I+2)-F1(I))/(X1(I+2)-X1(I))
4       XPI1=XP-X1(I+1)
        XPI=XP-X1(I)
        L1=XPI1*LP1
        L2=XPI*LP2
        P=F1(I)*(1.-2.*LP1*XPI)*L1*L1+F1(I+1)*(1.-2.*LP2*XPI1)
     *  *L2*L2+FP2*XPI1*L2*L2+FP1*XPI*L1*L1
        RETURN
        END
C-----------------------------------------------------------------------
C         YOU WILL NOT NEED THIS SUBROUTINE IF OTHER COMPILER
        SUBROUTINE GRAF(NP)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION IX(20000),IY(20000)
        COMMON /A/X1(20000),F1(20000)
        COMMON /B/SCAL1,SCAL2,IDL1,IDL2,XMIN,YMIN,XMAX,YMAX
        IDL1=0
        IDL2=430
        IPH1=639
        IPH2=60
        SCAL1=DBLE(IPH1-IDL1)/(XMAX-XMIN)
        SCAL2=DBLE(IDL2-IPH2)/(YMAX-YMIN)
        CALL VGA@
        CALL CLEAR_SCREEN@
        CALL SET_CURSOR_POS@(1,1)
        WRITE(5,'(A)')'ENTER = TO GO ON'
        CALL SET_CURSOR_POS@(40,1)
        WRITE(5,'(A14,E10.3)')'MAX. Y-VALUE =',YMAX
        CALL SET_CURSOR_POS@(40,28)
        WRITE(5,'(A14,E10.3)')'MIN. Y-VALUE =',YMIN
C               OS X
        CALL DRAW_LINE@(IDL1,IDL2,IPH1,IDL2,15)
        CALL DRAW_LINE@(IDL1,IPH2,IPH1,IPH2,15)
C               OS Y
        CALL DRAW_LINE@(IDL1,IDL2,IDL1,IPH2,15)
        CALL DRAW_LINE@(IPH1,IDL2,IPH1,IPH2,15)
        DO 20 J=1,NP
          IX(J)=DINT((X1(J)-XMIN)*SCAL1)+IDL1
          IY(J)=-DINT((F1(J)-YMIN)*SCAL2)+IDL2
20      CONTINUE
C               SPEKTRUM
        CALL POLYLINE@(IX,IY,NP,3)
        RETURN
        END
C-----------------------------------------------------------------------
C         YOU WILL NOT NEED THIS SUBROUTINE IF OTHER COMPILER
        SUBROUTINE GRAF1(XX1,YY1)
        COMMON /B/SCAL1,SCAL2,IDL1,IDL2,XMIN,YMIN,XMAX,YMAX
          IRX=DINT((XX1-XMIN)*SCAL1)+IDL1
          IRY=-DINT((YY1-YMIN)*SCAL2)+IDL2
          CALL SET_PIXEL@(IRX,IRY,14)
        RETURN
        END
C-----------------------------------------------------------------------
       SUBROUTINE HUNT(N,X,JLO)
C       TO FIND THE CLOSEST SMALER DATA POINT (JLO-TH) TO THE X-VALUE
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /A/ XX(20000),F1(20000)
C      DIMENSION XX(N)
      LOGICAL ASCND
      ASCND=XX(N).GT.XX(1)
      IF(JLO.LE.0.OR.JLO.GT.N)THEN
        JLO=0
        JHI=N+1
        GO TO 3
      ENDIF
      INC=1
      IF(X.GE.XX(JLO).EQV.ASCND)THEN
1       JHI=JLO+INC
        IF(JHI.GT.N)THEN
          JHI=N+1
        ELSE IF(X.GE.XX(JHI).EQV.ASCND)THEN
          JLO=JHI
          INC=INC+INC
          GO TO 1
        ENDIF
      ELSE
        JHI=JLO
2       JLO=JHI-INC
        IF(JLO.LT.1)THEN
          JLO=0
        ELSE IF(X.LT.XX(JLO).EQV.ASCND)THEN
          JHI=JLO
          INC=INC+INC
          GO TO 2
        ENDIF
      ENDIF
3     IF(JHI-JLO.EQ.1)RETURN
      JM=(JHI+JLO)/2
      IF(X.GT.XX(JM).EQV.ASCND)THEN
        JLO=JM
      ELSE
        JHI=JM
      ENDIF
      GO TO 3
      END
C-----------------------------------------------------------------------
        SUBROUTINE HUNT1(N,X,JLO)
C       TO FIND THE CLOSEST SMALER DATA POINT (JLO-TH) TO THE X-VALUE
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON /BET/ XX(1000),FLTRY(1000),NDAT2
C      DIMENSION XX(N)
      LOGICAL ASCND
      ASCND=XX(N).GT.XX(1)
      IF(JLO.LE.0.OR.JLO.GT.N)THEN
        JLO=0
        JHI=N+1
        GO TO 3
      ENDIF
      INC=1
      IF(X.GE.XX(JLO).EQV.ASCND)THEN
1       JHI=JLO+INC
        IF(JHI.GT.N)THEN
          JHI=N+1
        ELSE IF(X.GE.XX(JHI).EQV.ASCND)THEN
          JLO=JHI
          INC=INC+INC
          GO TO 1
        ENDIF
      ELSE
        JHI=JLO
2       JLO=JHI-INC
        IF(JLO.LT.1)THEN
          JLO=0
        ELSE IF(X.LT.XX(JLO).EQV.ASCND)THEN
          JHI=JLO
          INC=INC+INC
          GO TO 2
        ENDIF
      ENDIF
3     IF(JHI-JLO.EQ.1)RETURN
      JM=(JHI+JLO)/2
      IF(X.GT.XX(JM).EQV.ASCND)THEN
        JLO=JM
      ELSE
        JHI=JM
      ENDIF
      GO TO 3
      END