C **********************************************************************
C "EL2000C.F" - PROGRAM COMPUTING AN ORBIT OF A METEOR, WHICH GEOCENTRIC
C               RADIANT AND VELOCITY ARE KNOWN
C               VERSION: January 20, 2009
C             - TO EQUINOX 2000.0
C **********************************************************************
C
C MAIN PROGRAM
C
C   INPUT: all data on meteors in the standardized format
C
C   OUTPUT:
C      EA(1) - PERIHELION DISTANCE [AU]
C      EA(2) - ECCENTRICITY
C      EA(3) - ARGUMENT OF PERIHELION [DEGREES]
C      EA(4) - LONGITUDE OF NODE [DEGREES]
C      EA(5) - INCLINATION TO ECLIPTIC [DEGREES]
C      EA(6) - (empty)
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,K-Z)
      IMPLICIT INTEGER(I,J)
      CHARACTER*4 MON
      character*24 filin,fil15,fil16
      character*72 ret72
      DIMENSION EA(6),MON(12),EB(6),R(6),WEL(7),iau(5000),isol(5000)
      dimension sollam(5000),al(5000),dl(5000),vg(5000)
      dimension q(5000),e(5000),som(5000),gom(5000),sk(5000),jnn(5000)

      PI=4.0D0*DATAN(1.0D0)
      PI180=PI/0.180D3
      GAUSS=1.7202098950D-2
      G2=GAUSS*GAUSS
      GMS=0.2959122082855911E-03
      GMB=0.8997011346712499E-09
      AU=0.1495978706910000E+09  
      AUKM=0.149597870D9
      D24=8.6400D4
      VELU=D24/AUKM
C      EPS=2.344579D1*PI180 ... 1950.0
      EPS=2.3439291D1*PI180
      CEPS=DCOS(EPS)
      SEPS=DSIN(EPS)
C
      MON(1)='JAN.'
      MON(2)='FEB.'
      MON(3)='MAR.'
      MON(4)='APR.'
      MON(5)='MAY '
      MON(6)='JUNE'
      MON(7)='JULY'
      MON(8)='AUG.'
      MON(9)='SEP.'
      MON(10)='OCT.'
      MON(11)='NOV.'
      MON(12)='DEC.'
C
c Tolerance limits:
      open(unit=9,file='inparams.ele',access='sequential')
        read(9,*) ret72
        read(9,*) filin
        read(9,*) ret72
        read(9,*) fil15
        read(9,*) ret72
        read(9,*) fil16
        read(9,*) ret72
        read(9,*) DELQ
        read(9,*) ret72
        read(9,*) DELE
        read(9,*) ret72
        read(9,*) DELSO
        read(9,*) ret72
        read(9,*) DELGO
        read(9,*) ret72
        read(9,*) DELSK
      close(unit=9)

      open(unit=10,file=filin,access='sequential')
      j=0
  10  continue
      j=j+1
      read(10,*,end=20) iau(j),isol(j),sollam(j),al(j),dl(j),vg(j),q(j),
     *e(j),som(j),gom(j),sk(j),jnn(j)
      goto 10
  20  continue
      jall=j-1
      close(unit=10)

      OPEN(UNIT=15,FILE=fil15,ACCESS='SEQUENTIAL')
      OPEN(UNIT=16,FILE=fil16,ACCESS='SEQUENTIAL')
      WRITE(15,500)
      WRITE(15,510)
C
      DO 100 i=1,jall
        zemlam=sollam(i)+1.80d2
        if(zemlam.ge.3.60d2) zemlam=zemlam-3.60d2
        ttv=(zemlam*pi180-1.75347031435d0)/6.28307584918d3
        t20=ttv*1.0d1
        call EARTHrv(t20,X,Y,Z,VXE,VYE,VZE,WEL)
c        ET=3.65250d5*ttv+2.4515450d6
c        J=3
c        CALL  Pleph ( ET, J, NCTR, R )
C
cC The rectangular ecliptical coordinates of the Earth/meteor:
c        X=R(1)
c        Y=R(2)*CEPS+R(3)*SEPS
c        Z=-R(2)*SEPS+R(3)*CEPS
        RE=DSQRT(X*X+Y*Y+Z*Z)
cC The rectangular ecliptical components of the Earth's velocity vector:
c        VXE=R(4)
c        VYE=R(5)*CEPS+R(6)*SEPS
c        VZE=-R(5)*SEPS+R(6)*CEPS
cC
        VXR=VXE
        VYR=VYE*CEPS-VZE*SEPS
        VZR=VYE*SEPS+VZE*CEPS
        SDZ=VZR
        CDZ=DSQRT(1.0D0-SDZ*SDZ)
        ALZ=DATAN(VYR/VXR)
        ALZ=DABS(ALZ)
        IF(VXR.LT.0.0D0.AND.VYR.GE.0.0D0) ALZ=PI-ALZ
        IF(VXR.LT.0.0D0.AND.VYR.LT.0.0D0) ALZ=PI+ALZ
        IF(VXR.GE.0.0D0.AND.VYR.LT.0.0D0) ALZ=2.0D0*PI-ALZ
C
        YMIN=9.999D99
        YMAX=0.0D0
C
        CALA=DCOS(al(i)*PI180)
        SALA=DSIN(al(i)*PI180)
        CDLA=DCOS(dl(i)*PI180)
        SDLA=DSIN(dl(i)*PI180)
        VXG=vg(i)*VELU*CALA*CDLA
        VYR=vg(i)*VELU*SALA*CDLA
        VZR=vg(i)*VELU*SDLA
        VYG=VYR*CEPS+VZR*SEPS
        VZG=-VYR*SEPS+VZR*CEPS
        VX=-(VXG-VXE)
        VY=-(VYG-VYE)
        VZ=-(VZG-VZE)
C
        VVH=DSQRT(VX*VX+VY*VY+VZ*VZ)/VELU
C
        CALL ELEMS(t20,X,Y,Z,VX,VY,VZ,EB)
C
CC      CGAM=SDZ*SDLA+CDZ*CDLA*DCOS(ALZ-YY(I,11)*PI180)
CC      VE2L=8.8804D2*(2.0D0/RE-1.0D0)
CC      VEL=DSQRT(VE2L)
CC      CGAM=(VE2L+YY(I,14)*YY(I,14)-YY(I,15)*YY(I,15))/2.0D0/VEL/YY(I,14)
CC      SGAM=DSQRT(1.0D0-CGAM*CGAM)
CC      GAMMA=DATAN(SGAM/CGAM)
CC      IF(GAMMA.LT.0.0D0) GAMMA=PI+GAMMA
CC      GAMMA=GAMMA/PI180
C
 180  CONTINUE
        QW=EB(1)
        AW1=(1.0D0-EB(2))/EB(1)
        AW=1.0D0/AW1
        QQW=AW*(1.0D0+EB(2))
        EW=EB(2)
        SKW=EB(5)
        SOW=EB(3)
        GOW=EB(4)
        PIW=EB(3)+EB(4)
        IF(PIW.GE.3.60D2) PIW=PIW-3.60D2
CC      GAMW=GAMMA
        VINFW=DSQRT(vg(i)*vg(i)+1.249924D2)
        XX=VINFW
        IF(XX.LT.YMIN) YMIN=XX
        IF(XX.GT.YMAX) YMAX=XX
C      QS=QW
C      AS=AW
C      ES=EW
C      SKS=SKW
C      SOS=SOW
C      GOS=GOW
C      PIS=PIW
C      LAMS=GAMW
C      MNS=MNO
C      VVHS=VVH
C      VIS=VINFW
C  22  CONTINUE
C  44  CONTINUE
C  66  CONTINUE
      WRITE(15,610) i,iau(i),isol(i),q(i),e(i),som(i),gom(i),sk(i)
      WRITE(15,650) QW,EW,SOW,GOW,SKW
C
CC      IF(FQ(I,6).EQ.' 1') THEN
CC      IF(FE(I,6).EQ.' 1') THEN
CC      DELX=DYY(I,6)
CC      ELSE
CC      DELX=1.0D-2*YY(I,6)
CC      END IF
CC      IF(DABS(YY(I,6)-GAMW).GT.DELX) WRITE(16,*) CID(I),' [',ANO(I),']: 
CC     *  LS_obs,   LS_calc = ',YY(I,6),GAMW
CCCC      ELSE
CCCC      WRITE(16,*) CID(I),' [',ANO(I),']:  impossible checking of q; obse
CCCC     *rved value is missing .'
CC      END IF
      IF(DABS(q(i)-QW).GT.DELQ) WRITE(16,520) iau(i),isol(i),q(i),QW
      IF(DABS(e(i)-EW).GT.DELE) WRITE(16,530) iau(i),isol(i),e(i),EW
      IF(DABS(sk(i)-SKW).GT.DELSK) WRITE(16,540) iau(i),isol(i),sk(i),SK
     *W
      IF(DABS(som(i)-SOW).GT.DELSO) WRITE(16,550) iau(i),isol(i),som(i),
     *SOW
      IF(DABS(gom(i)-GOW).GT.DELGO) WRITE(16,560) iau(i),isol(i),gom(i),
     *GOW
 100  CONTINUE
      CLOSE(UNIT=10)
      CLOSE(UNIT=15)
      CLOSE(UNIT=16)
      write(*,*) 'Program successfully terminated.'
C
      STOP
C
 500  FORMAT('ser.No. IAU AdNo   q        e       arg.     node       i'
     *)
 510  FORMAT('----------------------------------------------------------
     *-')
 520  format('IAU No. =',i5,' AdNo. =',i3,'    q_orig =',f6.3,'    q_rec
     *. =',f6.3)
 530  format('IAU No. =',i5,' AdNo. =',i3,'    e_orig =',f6.3,'    e_rec
     *. =',f6.3)
 540  format('IAU No. =',i5,' AdNo. =',i3,'    i_orig =',f6.1,'    i_rec
     *. =',f6.1)
 550  format('IAU No. =',i5,' AdNo. =',i3,'  arg_orig =',f6.1,'  arg_rec
     *. =',f6.1)
 560  format('IAU No. =',i5,' AdNo. =',i3,'  nod_orig =',f6.1,'  nod_rec
     *. =',f6.1)
 600  FORMAT(A8,1X,A6,2X,I4,I3,2F9.5,F8.4,F9.5,2F9.3,F11.5)
 610  FORMAT(i5,i6,i3,2f9.5,3f9.3)
 620  FORMAT(A8,1X,A6,2X,I4,I3,2F9.5,8X,F9.5,2F9.3,F11.5)
 630  FORMAT(A8,1X,A6,2X,I4,I3,F9.5,9X,F8.4,F9.5,2F9.3,F11.5)
 632  FORMAT(A8,1X,A6,18X,F9.5,F8.4,F9.5,2F9.3,F11.5,/)
 650  FORMAT(14x,2f9.5,3f9.3,/)
C
      END 
C ---------------------------------------------------------------------
      SUBROUTINE ELEMS(T20,X,Y,Z,VX,VY,VZ,EA)
C ---------------------------------------------------------------------
C      - CALCULATION OF ORBIT OF A BODY, WHICH POSITION IN A GIVEN TIME
C        IS CHARACTERIZED WITH RADIUS AND VELOCITY VECTORS
C
C   INPUT QUANTITIES:
C        T20 - TIME, WHEN BODY'S POSITION IS KNOWN [JULIAN CENTURIES
C              FROM 2000.0; JD = 2451545.0]
C        X, Y, Z - RECTANGULAR ECLIPTIC HELIOCENTRIC COORDINATES OF
C                  THE BODY IN TIME T20 [AU]
C        VX, VY, VZ - COMPONENTS OF BODY'S VELOCITY VECTOR IN THE EC-
C                     LIPTIC HELIOCENTRIC COORDINATE SYSTEM IN TIME T20
C                     [AU/DAY]
C
C   OUTPUT QUAMTITIES (ELEMENTS OF BODY'S ORBIT):
C        EA(1) - PERIHELION DISTANCE [AU]
C        EA(2) - ECCENTRICITY
C        EA(3) - ARGUMENT OF PERIHELION [DEGREES]
C        EA(4) - LONGITUDE OF NODE [DEGREES]
C        EA(5) - INCLINATION TO ECLIPTIC [DEGREES]
C        EA(6) - TIME OF PERIHELION PASSAGE [JULIAN CENTURIES FROM
C                2000.0; JD = 2451545.0]
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EA(6)
C
      PI=4.0D0*DATAN(1.0D0)
      PI05=PI/2.0D0
      PI2=2.0D0*PI
      PI180=PI/0.180D3
      GAUSS=0.17202098950D-1
      G2=GAUSS*GAUSS
C
      R=DSQRT(X*X+Y*Y+Z*Z)
      V=DSQRT(VX*VX+VY*VY+VZ*VZ)
      CVFI=-(X*VX+Y*VY+Z*VZ)/R/V
      SVFI2=1.0D0-CVFI*CVFI
      YMEN=2.0D0/R-V*V/G2
      IF(DABS(YMEN).GT.1.0D-12) GOTO 10
      EA(1)=R*SVFI2
      EA(2)=1.0D0
      GOTO 20
  10  CONTINUE
      A=1.0D0/YMEN
      EA(2)=DSQRT(1.0D0-R*(2.0D0-R/A)*SVFI2/A)
      EA(1)=A*(1.0D0-EA(2))
C
  20  CONTINUE
      IF(EA(2).LT.1.0D-12) GOTO 40
      CF0=(EA(1)*(1.0D0+EA(2))/R-1.0D0)/EA(2)
      IF(DABS(CF0).GT.1.000001D0) WRITE(5,*) 'ERROR IN ELEMS; CF0'
      IF(DABS(CF0).GT.1.000001D0) STOP
      IF(CF0.GT.1.0D0) CF0=0.999999999D0
      IF(CF0.LT.-1.0D0) CF0=-0.999999999D0
      SF0=DSQRT(1.0D0-CF0*CF0)
      IF(DABS(CF0).LT.1.0D-12) GOTO 30
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=PI+F0
      F=F0
      IF(CVFI.GT.0.0D0) F=PI2-F0
      GOTO 40
  30  F=PI/2.0D0
C
  40  CONTINUE
      BX=Y*VZ-Z*VY
      BY=Z*VX-X*VZ
      BZ=X*VY-Y*VX
      CI=BZ/DSQRT(BX*BX+BY*BY+BZ*BZ)
      IF(DABS(CI).GT.1.000001D0) WRITE(5,*) ' ERROR IN SUBROUTINE ELEMS 
     *(CALCULATION OF "CI")'
      IF(DABS(CI).GT.1.000001D0) STOP
      IF(CI.GT.1.0D0) CI=0.999999999D0
      IF(CI.LT.-1.0D0) CI=-0.999999999D0
      SI=DSQRT(1.0D0-CI*CI)
      IF(DABS(CI).LT.1.0D-12) GOTO 42
      YI=DATAN(SI/CI)
      IF(YI.LT.0.0D0) YI=PI+YI
      EA(5)=YI/PI180
      GOTO 45
  42  EA(5)=9.0D1
C
  45  CONTINUE
      CGO=-BY/DSQRT(BX*BX+BY*BY)
      SGO=BX/DSQRT(BX*BX+BY*BY)
      IF(DABS(CGO).LT.1.0D-12) GOTO 47
      YGO=DATAN(SGO/CGO)
      YGO=DABS(YGO)
      GOTO 50
  47  YGO=PI/2.0D0
  50  CONTINUE
      IF(CGO.GE.0.0D0.AND.SGO.LT.0.0D0) YGO=PI2-YGO
      IF(CGO.LT.0.0D0.AND.SGO.GE.0.0D0) YGO=PI-YGO
      IF(CGO.LT.0.0D0.AND.SGO.LT.0.0D0) YGO=PI+YGO
      EA(4)=YGO/PI180
C
      IF(EA(2).LT.1.0D-12) GOTO 60
      CSOF=(X*CGO+Y*SGO)/R
      IF(DABS(CI).LT.1.0D-12) GOTO 52
      SSOF=(Y*CGO-X*SGO)/R/CI
      GOTO 54
  52  SSOF=Z/R/SI
  54  CONTINUE
      IF(DABS(CSOF).LT.1.0D-12) GOTO 56
      SOF=DATAN(SSOF/CSOF)
      SOF=DABS(SOF)
      GOTO 58
  56  SOF=PI/2.0D0
  58  CONTINUE
      IF(CSOF.GE.0.0D0.AND.SSOF.LT.0.0D0) SOF=PI2-SOF
      IF(CSOF.LT.0.0D0.AND.SSOF.GE.0.0D0) SOF=PI-SOF
      IF(CSOF.LT.0.0D0.AND.SSOF.LT.0.0D0) SOF=PI+SOF
      YSO=SOF-F
      IF(YSO.LT.0.0D0) YSO=YSO+PI2
      EA(3)=YSO/PI180
      GOTO 70
  60  CONTINUE
      EA(3)=0.0D0
C
  70  CONTINUE
C
      RETURN
C
      END
C ----------------------------------------------------------------------
      SUBROUTINE EARTHrv(T20,X,Y,Z,VXZ,VYZ,VZZ,WEL)
C ----------------------------------------------------------------------
C     - CALCULATION OF QUANTITIES CHARACTERIZING THE POSITION OF THE
C       EARTH AND ITS ORBIT AT GIVEN TIME (RECTANGULAR COORDINATES OF
C       THE EARTH WERE CALCULATED BY: BRETAGNON P.: 1982, ASTRON.
C       ASTROPHYS. 114, 278.)
C
C INPUT:
C    T20 - THE TIME [JULIAN CENTURIES FROM JD = 2451545.0]
C
C OUTPUT = THE VALUES VALID IN TIME "T20" (EQUINOX 2000.0):
C    X, Y, Z - HELIOCENTRIC RECTANGULAR COORDINATES REFERRED TO THE
C              ECLIPTIC [AU]
C    WEL(1) - SEMI-MAJOR AXIS [AU]
C    WEL(2) - ECCENTRICITY
C    WEL(3) - ECLIPTIC MEAN LONGITUDE [RADIANS]
C    WEL(4) - LONGITUDE OF PERIHELION [RADIANS]
C    WEL(5) - LONGITUDE OF ASCENDING NODE [RADIANS]
C    WEL(6) - INCLINATION TO THE ECLIPTIC [RADIANS]
C    WEL(7) - TRUE ANOMALY [RADIANS]
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8 (A-H,K-Z)
      DIMENSION WEL(7)
C
      PI=4.0D0*DATAN(1.0D0)
      PI2=2.0D0*PI
      GAUSS=1.7202098950D-2
C
      WEL(1)=1.00000101778D0
C
      WEL(3)=1.75347031435D0+6.28307584918D2*T20
      WEL(3)=WEL(3)-9.9189D-8*T20*T20+7.3D-13*T20*T20*T20
   3  CONTINUE
      IF(WEL(3).GE.0.0D0) GOTO 5
      WEL(3)=WEL(3)+PI2
      GOTO 3
   5  CONTINUE
      IF(WEL(3).LT.PI2) GOTO 7
      WEL(3)=WEL(3)-PI2
      GOTO 5
C
   7  WK=-3.7408165D-3-8.2266699D-5*T20
      WK=WK+2.748939D-7*T20*T20+1.04217D-9*T20*T20*T20
      WH=1.628447663D-2-6.2030259D-5*T20
      WH=WH-3.353888D-7*T20*T20+7.1185D-10*T20*T20*T20
      WQ=-1.13469002D-4*T20+1.237314D-7*T20*T20+1.2705D-9*T20*T20*T20
      WP=1.0180391D-5*T20+4.701998D-7*T20*T20-5.3829D-10*T20*T20*T20
C 
      WEL(2)=DSQRT(WK*WK+WH*WH)
C
      IF(DABS(WK).LT.1.0D-12) GOTO 210
      WEL(4)=DATAN(WH/WK)
      WEL(4)=DABS(WEL(4))
      IF(WH.GE.0.0D0.AND.WK.LT.0.0D0) WEL(4)=PI-WEL(4)
      IF(WH.LT.0.0D0.AND.WK.LT.0.0D0) WEL(4)=PI+WEL(4)
      IF(WH.LT.0.0D0.AND.WK.GE.0.0D0) WEL(4)=PI2-WEL(4)
C
  10  CONTINUE
      IF(DABS(WQ).LT.1.0D-12) GOTO 220
      WEL(5)=DATAN(WP/WQ)
      WEL(5)=DABS(WEL(5))
      IF(WP.GE.0.0D0.AND.WQ.LT.0.0D0) WEL(5)=PI-WEL(5)
      IF(WP.LT.0.0D0.AND.WQ.LT.0.0D0) WEL(5)=PI+WEL(5)
      IF(WP.LT.0.0D0.AND.WQ.GE.0.0D0) WEL(5)=PI2-WEL(5)
  20  WG=DSQRT(WP*WP+WQ*WQ)
      POM=1.0D0-WG*WG
      IF(DABS(POM).LT.1.0D-12) GOTO 230
      POM=WG/DSQRT(POM)
      WEL(6)=2.0D0*DATAN(POM)
      IF(WEL(6).LT.0.0D0) WEL(6)=PI+WEL(6)
C
  30  CONTINUE
      WM=WEL(3)-WEL(4)
      WEL(7)=WM+(2.0D0-WEL(2)*WEL(2)/4.0D0)*WEL(2)*DSIN(WM)
      WEL(7)=WEL(7)+5.0D0*WEL(2)*WEL(2)/4.0D0*DSIN(2.0D0*WM)
      WEL(7)=WEL(7)+1.3D1*WEL(2)*WEL(2)*WEL(2)/1.2D1*DSIN(3.0D0*WM)
  40  CONTINUE
      IF(WEL(7).GE.0.0D0) GOTO 50
      WEL(7)=WEL(7)+PI2
      GOTO 40
  50  CONTINUE
      IF(WEL(7).LT.PI2) GOTO 60
      WEL(7)=WEL(7)-PI2
      GOTO 50
  60  RZ=WEL(1)*(1.0D0-WEL(2)*WEL(2))/(1.0D0+WEL(2)*DCOS(WEL(7)))
      W2=WEL(7)+WEL(4)-WEL(5)
      C2O=DCOS(W2)
      S2O=DSIN(W2)
      CVO=DCOS(WEL(5))
      SVO=DSIN(WEL(5))
      CSK=DCOS(WEL(6))
      SSK=DSIN(WEL(6))
      X=RZ*(C2O*CVO-S2O*CSK*SVO)
      Y=RZ*(C2O*SVO+S2O*CSK*CVO)
      Z=RZ*S2O*SSK
      VELZ=GAUSS*DSQRT(2.0D0/RZ-1.0D0/WEL(1))

      MEN2=1.0D0+WEL(2)*WEL(2)+2.0D0*WEL(2)*DCOS(wel(7))
      MEN=DSQRT(MEN2)
      CVFI=-WEL(2)*DSIN(WEL(7))/MEN
      SVFI=(1.0D0+WEL(2)*DCOS(WEL(7)))/MEN
      VFI=DATAN(SVFI/CVFI)
      VFI=DABS(VFI)
      if(CVFI.LT.0.0D0.AND.SVFI.GE.0.0D0) VFI=PI-VFI
      if(CVFI.LT.0.0D0.AND.SVFI.LT.0.0D0) VFI=PI+VFI
      if(CVFI.GE.0.0D0.AND.SVFI.LT.0.0D0) VFI=PI2-VFI
      SOM=WEL(4)-WEL(5)
      C3Q=DCOS(SOM+WEL(7)-VFI)
      S3Q=DSIN(SOM+WEL(7)-VFI)
      CGO=DCOS(WEL(5))
      SGO=DSIN(WEL(5))
      CSK=DCOS(WEL(6))
      SSK=DSIN(WEL(6))
      VXZ=-VELZ*(C3Q*CGO-S3Q*CSK*SGO)
      VYZ=-VELZ*(C3Q*SGO+S3Q*CSK*CGO)
      VZZ=-VELZ*S3Q*SSK
      RETURN
C
 210  WEL(4)=PI/2.0D0
      IF(WH.LT.0.0D0) WEL(4)=1.5D0*PI
      GOTO 10
 220  WEL(5)=PI/2.0D0
      IF(WP.LT.0.0D0) WEL(5)=1.5D0*PI
      GOTO 20
 230  WEL(6)=PI
      GOTO 30
C
      END
