C ********************************************************************** C "METHODS.FOR" - PROGRAM COMPUTING THEORETICAL METEOR SHOWER RADIANTS C (VERSION 2) C - TO EQUINOX 2000.0 C ********************************************************************** C C STATUS OF USING: C ==================== C C 1. THE USE OF THE PROGRAM AS A SUBJECT OF COMMERCE IS NOT ALLOWED. C 2. FOR OTHER PURPOSES, THE AUTHORS PERMIT TO COPY AND USE THE C PROGRAM OR ITS PART(S) FREELY; C 3. EACH USING OF THE PROGRAM OR ANY OF ITS PARTS (SUBROUTINES) C SHOULD BE QUOTED IN THE OFFICIAL PRESENTATION; THE PROGRAM WAS C INTRODUCED IN THE "ASTRONOMY AND ASTROPHYSICS, LETTERS", PLEASE C SEE: NESLUSAN L., SVOREN J., PORUBCAN V.: ASTRON. ASTROPHYS., C SUBMITTED. C C ********************************************************************** C C BRIEF DESCRIPTION: C ================== C C THE PROGRAM COMPUTES A THEORETICAL RADIANT OF A METEOR SHOWER, C IF THE ORBIT OF AN ASSUMED PARENT BODY IS WELL-KNOWN C C THE PROGRAM INCLUDES SIX DIFFERENT METHODS: C - (Q) Q-ADJUSTMENT (HASEGAWA; 1990) C - (B) VARIATION OF PERIHELION DISTANCE AND EXCENTRICITY (SVOREN C ET AL; 1993) C - (W) ROTATION OF THE LINE OF APSIDES (STEEL AND BAGGALEY; 1985) C - (A) ROTATION AROUND THE LINE OF APSIDES (SVOREN ET AL; 1993) C - (H) OMEGA-ADJUSTMENT (HASEGAWA; 1990) C - (P) PORTER'S METHOD (PORTER; 1952) C (FOR MORE DETAILED DESCRIPTION AND REFERENCES, PLEASE, SEE IN THE C ABOVE MENTIONED PAPERS. REFERENCES ARE ALSO GIVEN IN BRIEF INDI- C VIDUAL DESCRIPTIONS OF SUBROUTINES BELOW.) C C --------------------------------------------------------------------- C MAIN PROGRAM C C INPUT: C ORBITAL ELEMENTS OF THE ORBIT OF THE PARENT BODY 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 YEAR, MONTH, DAY - TIME OF PERIHELION PASSAGE OF THE PARENT BODY C (CONVERTED TO: EA(6) - THE TIME IN JULIAN CENTURIES FROM C JD = 2433282.423) C C OUTPUT: C THE QUANTITIES DESCRIBING THEORETICAL RADIANTS AT BOTH ARCS (PRE- C PERIHELION AND POST-PERIHELION) OF THE ORBIT OF THE PARENT BODY; C THE OUTPUT QUANTITIES AT THE POST-PERIHELION ARC ARE: C R1(1) - SOLAR LONGITUDE OF THE MAXIMUM OF POTENTIAL SHOWER C [DEGRRES] C R1(2) - RIGHT ASCENSION OF THE PREDICTED RADIANT [DEGREES] C R1(3) - DECLINATION OF THE RADIANT [DEGREES] C R1(4) - PREDICTED GEOCENTRIC VELOCITY OF METEORS [KM/S] C R1(5) - PREDICTED HELIOCENTRIC VELOCITY OF METEORS [KM/S] C R1(6) - PREDICTED TIME OF THE SHOWER MAXIMUM [JULIAN CENTURIES C FROM 1950.0] C WM1(1) - MOMENT OF MINIMUM DISTANCE BETWEEN THE EARTH AND THE C ORBIT OF THE PARENT BODY [JULIAN CENTURIES FROM 1950.0] C WM1(2) - MINIMUM DISTANCE BETWEEN POST-PERIHELION ARC OF THE C PARENT BODY ORBIT AND EARTH'S ORBIT [AU] C WM1(3) - DISTANCE OF THE PARENT BODY FROM THE EARTH AT THE MOMENT C OF MAXIMUM [AU] C WM1(4) - TIME INTERVAL BETWEEN THE PASSAGES OF THE PARENT BODY C AND EARTH ACROSS THE NEAREST POINTS OF THE POST-PERI- C HELION ARC OF THE PARENT BODY ORBIT AND EARTH'S ORBIT C [DAYS] C WM1(5) - TRUE ANOMALY OF THE PARENT BODY AT THE MOMENT OF ITS C NEAREST POST-PERIHELION APPROACH TO THE EARTH'S ORBIT C [DEGREES] C WM1(6) - TRUE ANOMALY OF THE EARTH AT THE MOMENT OF ITS NEAREST C APPROACH TO THE POST-PERIHELION ARC OF THE PARENT BODY C ORBIT [DEGREES] C ORBITAL ELEMENTS OF MODIFIED ORBIT CROSSING THE EARTH'S ORBIT C EB1(1) - PERIHELION DISTANCE [AU] C EB1(2) - ECCENTRICITY C EB1(3) - ARGUMENT OF PERIHELION [DEGREES] C EB1(4) - LONGITUDE OF NODE [DEGREES] C EB1(5) - INCLINATION TO ECLIPTIC [DEGREES] C DD1 - D-CRITERION CHARACTERIZING THE DIFFERENCE BETWEEN THE ORBIT C OF THE PARENT BODY (INPUT ORBIT) AND MODIFIED ORBIT C CROSSING THE EARTH'S ORBIT (OUTPUT ORBIT) C THE ANALOGOUS QUANTITIES CONCERNING THE PRE-PERIHELION ARC ARE C STORED IN THE FIELDS OF VALUES "R2(I)", "WM2(I)", AND "EB2(I)", C WHERE I=1, 2,..., 6 (5, RESPECTIVELY), AND IN VALUE "DD2". C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) INTEGER I,IM,IYR,J C CHARACTER*1 MET character*4 MON character*11 prm character*24 filin,fil15,fil16,fil17 character*72 ret72 DIMENSION EA(6),MON(12),EB(6),MET(6) DIMENSION EB1(6),EB2(6),R1(6),R2(6),RY1(6,6),RY2(6,6) DIMENSION EY1(6,6),EY2(6,6),DDY1(6),DDY2(6) DIMENSION WM1(6),WM2(6),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) C pi=4.0d0*datan(1.0d0) pi180=pi/1.80d2 c Tolerance limits: open(unit=9,file='inparams.rad',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,*) DELLS read(9,*) ret72 read(9,*) DELAL read(9,*) ret72 read(9,*) DELDL read(9,*) ret72 read(9,*) DELVG read(9,*) ret72 read(9,*) iyn close(unit=9) 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 MET(1)='Q' MET(2)='B' MET(3)='W' MET(4)='A' MET(5)='H' MET(6)='P' 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') if(iyn.eq.1) OPEN(UNIT=17,FILE='debug.d',ACCESS='SEQUENTIAL') WRITE(15,500) WRITE(15,510) IROKI=2000 DEN=1.0D0 IMES=1 CALL JULIAN(IROKI,IMES,DEN,T19,TMIN) IROKI=IROKI+1 CALL JULIAN(IROKI,IMES,DEN,T19,TMAX) IROKI=IROKI-1 DO 100 i=1,jall WRITE(*,*) 'METEOR SHOWER - IAU No., AdNo. = ',iau(i),isol(i) if(iyn.eq.1) WRITE(17,700) iau(i),isol(i) EA(1)=q(i) EA(2)=e(i) EA(3)=som(i) EA(4)=gom(i) EA(5)=sk(i) C c IROK=YY(I,3) c IMES=YY(I,4) c DEN=YY(I,5) cClll DEN=17.0 c CALL JULIAN(IROK,IMES,DEN,T19,ETT) c EA(6)=ETT CALL APPROACH(TMIN,TMAX,EA,WM1,WM2) if(iyn.eq.1) WRITE(17,590) WM1(2),WM2(2) zemlam=sollam(i)+1.80d2 if(zemlam.ge.3.60d2) zemlam=zemlam-3.60d2 ttv=(zemlam*pi180-1.75347031435d0)/6.28307584918d3 ETT=3.65250d5*ttv+2.4515450d6 210 CONTINUE CALL QMETH(ETT,EA,R1,EB1,DD1,R2,EB2,DD2) Ct WRITE(*,*) 'R1(6) = ',R1(6) Ct WRITE(*,*) 'R1(1), GO = ',R1(1),EA(4) Ct WRITE(*,*) 'Von z QMETH.' if(iyn.eq.1) then WRITE(17,*) ' Q-method:' prm='D_1, D_2: ' WRITE(17,520) prm,DD1,DD2 prm='q: ' WRITE(17,530) prm,EA(1),EB1(1),EB2(1) prm='e: ' WRITE(17,530) prm,EA(2),EB1(2),EB2(2) prm='arg.: ' WRITE(17,530) prm,EA(3),EB1(3),EB2(3) prm='node: ' WRITE(17,530) prm,EA(4),EB1(4),EB2(4) prm='i: ' WRITE(17,530) prm,EA(5),EB1(5),EB2(5) prm='lambda_sun:' WRITE(17,530) prm,sollam(i),R1(1),R2(1) prm='R.A.: ' WRITE(17,530) prm,al(i),R1(2),R2(2) prm='DEC.: ' WRITE(17,530) prm,dl(i),R1(3),R2(3) prm='V_g: ' WRITE(17,530) prm,vg(i),R1(4),R2(4) prm='V_h(calc.):' WRITE(17,520) prm,R1(5),R2(5) endif DO 215 J=1,6 RY1(1,J)=R1(J) RY2(1,J)=R2(J) EY1(1,J)=EB1(J) EY2(1,J)=EB2(J) 215 CONTINUE DDY1(1)=DD1 DDY2(1)=DD2 Cmmm J1=1 Cmmm J2=1 Cmmm GOTO 283 C 220 CONTINUE CALL BMETH(ETT,EA,R1,EB1,DD1,R2,EB2,DD2) Ct WRITE(*,*) 'Von z BMETH.' if(iyn.eq.1) then WRITE(17,*) ' B-method:' prm='D_1, D_2: ' WRITE(17,520) prm,DD1,DD2 prm='q: ' WRITE(17,530) prm,EA(1),EB1(1),EB2(1) prm='e: ' WRITE(17,530) prm,EA(2),EB1(2),EB2(2) prm='arg.: ' WRITE(17,530) prm,EA(3),EB1(3),EB2(3) prm='node: ' WRITE(17,530) prm,EA(4),EB1(4),EB2(4) prm='i: ' WRITE(17,530) prm,EA(5),EB1(5),EB2(5) prm='lambda_sun:' WRITE(17,530) prm,sollam(i),R1(1),R2(1) prm='R.A.: ' WRITE(17,530) prm,al(i),R1(2),R2(2) prm='DEC.: ' WRITE(17,530) prm,dl(i),R1(3),R2(3) prm='V_g: ' WRITE(17,530) prm,vg(i),R1(4),R2(4) prm='V_h(calc.):' WRITE(17,520) prm,R1(5),R2(5) endif DO 225 J=1,6 RY1(2,J)=R1(J) RY2(2,J)=R2(J) EY1(2,J)=EB1(J) EY2(2,J)=EB2(J) 225 CONTINUE DDY1(2)=DD1 DDY2(2)=DD2 C 230 CONTINUE CALL WMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2) Ct WRITE(*,*) 'Von z WMETH.' if(iyn.eq.1) then WRITE(17,*) ' W-method:' prm='D_1, D_2: ' WRITE(17,520) prm,DD1,DD2 prm='q: ' WRITE(17,530) prm,EA(1),EB1(1),EB2(1) prm='e: ' WRITE(17,530) prm,EA(2),EB1(2),EB2(2) prm='arg.: ' WRITE(17,530) prm,EA(3),EB1(3),EB2(3) prm='node: ' WRITE(17,530) prm,EA(4),EB1(4),EB2(4) prm='i: ' WRITE(17,530) prm,EA(5),EB1(5),EB2(5) prm='lambda_sun:' WRITE(17,530) prm,sollam(i),R1(1),R2(1) prm='R.A.: ' WRITE(17,530) prm,al(i),R1(2),R2(2) prm='DEC.: ' WRITE(17,530) prm,dl(i),R1(3),R2(3) prm='V_g: ' WRITE(17,530) prm,vg(i),R1(4),R2(4) prm='V_h(calc.):' WRITE(17,520) prm,R1(5),R2(5) endif DO 235 J=1,6 RY1(3,J)=R1(J) RY2(3,J)=R2(J) EY1(3,J)=EB1(J) EY2(3,J)=EB2(J) 235 CONTINUE DDY1(3)=DD1 DDY2(3)=DD2 C 240 CONTINUE CALL AMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2) Ct WRITE(*,*) 'Von z AMETH.' if(iyn.eq.1) then WRITE(17,*) ' A-method:' prm='D_1, D_2: ' WRITE(17,520) prm,DD1,DD2 prm='q: ' WRITE(17,530) prm,EA(1),EB1(1),EB2(1) prm='e: ' WRITE(17,530) prm,EA(2),EB1(2),EB2(2) prm='arg.: ' WRITE(17,530) prm,EA(3),EB1(3),EB2(3) prm='node: ' WRITE(17,530) prm,EA(4),EB1(4),EB2(4) prm='i: ' WRITE(17,530) prm,EA(5),EB1(5),EB2(5) prm='lambda_sun:' WRITE(17,530) prm,sollam(i),R1(1),R2(1) prm='R.A.: ' WRITE(17,530) prm,al(i),R1(2),R2(2) prm='DEC.: ' WRITE(17,530) prm,dl(i),R1(3),R2(3) prm='V_g: ' WRITE(17,530) prm,vg(i),R1(4),R2(4) prm='V_h(calc.):' WRITE(17,520) prm,R1(5),R2(5) endif DO 245 J=1,6 RY1(4,J)=R1(J) RY2(4,J)=R2(J) EY1(4,J)=EB1(J) EY2(4,J)=EB2(J) 245 CONTINUE DDY1(4)=DD1 DDY2(4)=DD2 C 250 CONTINUE CALL HMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2) Ct WRITE(*,*) 'Von z HMETH.' if(iyn.eq.1) then WRITE(17,*) ' H-method:' prm='D_1, D_2: ' WRITE(17,520) prm,DD1,DD2 prm='q: ' WRITE(17,530) prm,EA(1),EB1(1),EB2(1) prm='e: ' WRITE(17,530) prm,EA(2),EB1(2),EB2(2) prm='arg.: ' WRITE(17,530) prm,EA(3),EB1(3),EB2(3) prm='node: ' WRITE(17,530) prm,EA(4),EB1(4),EB2(4) prm='i: ' WRITE(17,530) prm,EA(5),EB1(5),EB2(5) prm='lambda_sun:' WRITE(17,530) prm,sollam(i),R1(1),R2(1) prm='R.A.: ' WRITE(17,530) prm,al(i),R1(2),R2(2) prm='DEC.: ' WRITE(17,530) prm,dl(i),R1(3),R2(3) prm='V_g: ' WRITE(17,530) prm,vg(i),R1(4),R2(4) prm='V_h(calc.):' WRITE(17,520) prm,R1(5),R2(5) endif DO 255 J=1,6 RY1(5,J)=R1(J) RY2(5,J)=R2(J) EY1(5,J)=EB1(J) EY2(5,J)=EB2(J) 255 CONTINUE DDY1(5)=DD1 DDY2(5)=DD2 C 260 CONTINUE CALL PMETH(WM1,WM2,EA,R1,EB1,DD1,R2,EB2,DD2) Ct WRITE(*,*) 'Von z PMETH.' if(iyn.eq.1) then WRITE(17,*) ' P-method:' prm='D_1, D_2: ' WRITE(17,520) prm,DD1,DD2 prm='q: ' WRITE(17,530) prm,EA(1),EB1(1),EB2(1) prm='e: ' WRITE(17,530) prm,EA(2),EB1(2),EB2(2) prm='arg.: ' WRITE(17,530) prm,EA(3),EB1(3),EB2(3) prm='node: ' WRITE(17,530) prm,EA(4),EB1(4),EB2(4) prm='i: ' WRITE(17,530) prm,EA(5),EB1(5),EB2(5) prm='lambda_sun:' WRITE(17,530) prm,sollam(i),R1(1),R2(1) prm='R.A.: ' WRITE(17,530) prm,al(i),R1(2),R2(2) prm='DEC.: ' WRITE(17,530) prm,dl(i),R1(3),R2(3) prm='V_g: ' WRITE(17,530) prm,vg(i),R1(4),R2(4) prm='V_h(calc.):' WRITE(17,520) prm,R1(5),R2(5) WRITE(17,*) '' endif DO 265 J=1,6 RY1(6,J)=R1(J) RY2(6,J)=R2(J) EY1(6,J)=EB1(J) EY2(6,J)=EB2(J) 265 CONTINUE DDY1(6)=DD1 DDY2(6)=DD2 C DD1=DDY1(1) J1=1 DO 270 J=2,6 IF(DDY1(J).LT.DD1) J1=J IF(DDY1(J).LT.DD1) DD1=DDY1(J) 270 CONTINUE DD2=DDY2(1) J2=1 DO 280 J=2,6 IF(DDY2(J).LT.DD2) J2=J IF(DDY2(J).LT.DD2) DD2=DDY2(J) 280 CONTINUE C 283 CONTINUE IF(DD1.LT.9.9999D1.OR.DD2.LT.9.9999D1) GOTO 80 YQ=0.0D0 YE=0.0D0 YSO=0.0D0 YGO=0.0D0 YI=0.0D0 YLAM=0.0D0 ALB=0.0D0 DLB=0.0D0 VGB=0.0D0 VHB=0.0D0 GOTO 180 C 80 CONTINUE C ??? IF(DABS(RY2(J2,6)-EA(6)).LT.DABS(RY1(J1,6)-EA(6))) GOTO 135 IF(DD2.LT.DD1) GOTO 135 TJD=RY1(J1,6) CALL INVJUL(TJD,IRK,IMS,DNI) Ct WRITE(*,*) 'TJD, IRK, IMS, DNI = ',TJD,IRK,IMS,DNI YQ=EY1(J1,1) YE=EY1(J1,2) YSO=EY1(J1,3) YGO=EY1(J1,4) YI=EY1(J1,5) YLAMB=RY1(J1,1) ALB=RY1(J1,2) DLB=RY1(J1,3) VGB=RY1(J1,4) VHB=RY1(J1,5) JMET=J1 GOTO 180 135 CONTINUE TJD=RY2(J2,6) CALL INVJUL(TJD,IRK,IMS,DNI) YQ=EY2(J2,1) YE=EY2(J2,2) YSO=EY2(J2,3) YGO=EY2(J2,4) YI=EY2(J2,5) YLAMB=RY2(J2,1) ALB=RY2(J2,2) DLB=RY2(J2,3) VGB=RY2(J2,4) VHB=RY2(J2,5) JMET=J2 C 180 CONTINUE YLAMA=gom(i) IF(DABS(YLAMA-YLAMB).GT.1.80D2) YLAMA=YLAMA+1.80D2 IF(YLAMA.GE.3.60D2) YLAMA=YLAMA-3.60D2 WRITE(15,600) i,iau(i),isol(i),YLAMA,al(i),dl(i),vg(i) WRITE(15,650) YLAMB,ALB,DLB,VGB,MET(JMET) C lamx=sollam(i) if(lamx.lt.9.0d1.and.YLAMB.gt.2.70d2) lamx=lamx+3.60d2 if(lamx.gt.2.70d2.and.YLAMB.lt.9.0d1) lamx=lamx-3.60d2 IF(DABS(YLAMB-lamx).GT.DELLS) THEN WRITE(16,540) iau(i),isol(i),sollam(i),YLAMB END IF alx=al(i) if(alx.lt.9.0d1.and.ALB.gt.2.70d2) alx=alx+3.60d2 if(alx.gt.2.70d2.and.ALB.lt.9.0d1) alx=alx-3.60d2 IF(DABS(ALB-alx).GT.DELAL) THEN WRITE(16,550) iau(i),isol(i),al(i),ALB END IF IF(DABS(DLB-dl(i)).GT.DELDL) THEN WRITE(16,560) iau(i),isol(i),dl(i),DLB END IF IF(DABS(VGB-vg(i)).GT.DELVG) THEN WRITE(16,570) iau(i),isol(i),vg(i),VGB END IF c IF(DABS(VHB-vh(i)).GT.DELVG) THEN c WRITE(16,*) iau(i),isol(i),' Vh_obs, Vh_calc = ',vh(i),VHB c END IF CC CLOSE(UNIT=15) C 95 CONTINUE 100 CONTINUE CLOSE(UNIT=15) CLOSE(UNIT=16) if(iyn.eq.1) CLOSE(UNIT=17) write(*,*) 'Program successfully terminated.' C STOP C 500 FORMAT('ser.No. IAU AdNo LS RA DEC Vg method') 510 FORMAT('-------------------------------------------------------') 520 format(a11,8x,2f8.3) 530 format(a11,3f8.3) 540 format('IAU No. =',i5,' AdNo. =',i3,' LS_orig =',f7.2,' LS_re *c. =',f7.2) 550 format('IAU No. =',i5,' AdNo. =',i3,' RA_orig =',f7.2,' RA_re *c. =',f7.2) 560 format('IAU No. =',i5,' AdNo. =',i3,' DEC_orig =',f7.2,' DEC_re *c. =',f7.2) 570 format('IAU No. =',i5,' AdNo. =',i3,' V_g_orig =',f7.3,' V_g_re *c. =',f7.3) 590 format(' MOID: ',8x,2f8.4) 600 FORMAT(i5,i6,i3,3f9.3,f8.3) 610 FORMAT(A8,1X,A6,2X,I4,I3,F9.5,3F9.3,F8.3) 620 FORMAT(A8,1X,A6,2X,I4,I3,F9.5,9X,2F9.3,2F8.3) 650 FORMAT(14x,3f9.3,f8.3,3x,a1,/) 700 format('IAU No. =',i4,' AdNo. =',i3) END C --------------------------------------------------------------------- SUBROUTINE QMETH(ETT,EA,R1,EB1,DD1,R2,EB2,DD2) C --------------------------------------------------------------------- C - (Q) Q-ADJUSTMENT C (HASEGAWA I.: 1990, PUBL. ASTRON. SOC. JAPAN 42, 175.) C C "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE C POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 1950.0. NEXT C INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL- C DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2, C 3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS C THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D- C DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC C OF THE PARENT BODY ORBIT, RESPECTIVELY. C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) INTEGER I,IM,IYR,J DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7),EZ(6) C PI=4.0D0*DATAN(1.0D0) PI2=2.0D0*PI PI180=PI/0.18D3 EPS=2.3439291D1 CEPS=DCOS(PI180*EPS) SEPS=DSIN(PI180*EPS) NTARG=3 NCENT=11 C DO 50 J=2,5 EB1(J)=EA(J) EB2(J)=EA(J) 50 CONTINUE DT0=(TMAX-TMIN)/3.6525D2 DT05=DT0/2.0D0 DD1=0.999D99 DD2=0.999D99 C IF(EA(3).GT.0.18D3) GOTO 100 FB=0.18D3-EA(3) IF(FB.GT.1.799999999D2.AND.EA(2).GT.0.9999999D0) GOTO 200 VO0=EA(4)+0.18D3 IF(VO0.GE.0.36D3) VO0=VO0-0.36D3 GOTO 110 100 FB=0.36D3-EA(3) VO0=EA(4) 110 VO0=VO0*PI180 DWM=2.0D0*PI2 TA=TMIN TB=TMAX DT=DT0 C 120 CONTINUE C>>> DO 130 ETT=TA,TB,DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) c XZ=RRD(1) c YZ=RRD(2)*CEPS+RRD(3)*SEPS C ZZ=-RRD(2)*SEPS+RRD(3)*CEPS C VXZ=RRD(4) C VYZ=RRD(5)*CEPS+RRD(6)*SEPS C VZZ=-RRD(5)*SEPS+RRD(6)*CEPS C CALL ELEMS(ETT,XZ,YZ,VZ,VXZ,VYZ,VZZ,EZ) WL1=WEL(7)+WEL(4) IF(WL1.GE.PI2) WL1=WL1-PI2 WL1=(WL1+PI)/PI180 CC WRITE(*,*) 'Old value = ',WL1 CCC T15=(ETT-2.433282423D6)/3.6525D4 CCC 121 CONTINUE CCC IF(WEL(3).LT.PI2) GOTO 123 CCC WEL(3)=WEL(3)-PI2 CCC GOTO 121 CCC 123 CONTINUE C IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF Ct WLOUT=(WL+PI)/PI180 Ct IF(WLOUT.GE.3.60D2) WLOUT=WLOUT-3.60D2 Ct WRITE(*,*) 'JPL eph., lam = ',WLOUT CCC WLNEW=EZ(6)+EZ(3)+EZ(4) CCC IF(WLNEW.GE.PI2) WLNEW=WLNEW-PI2 CCC IF(WLNEW.GE.PI2) WLNEW=WLNEW-PI2 CCC WLNEW=(WLNEW+PI)/PI180 CCC WRITE(*,*) 'New value = ',WLNEW CCC WRITE(*,*) ' ' CCC STOP C>>> IF(DABS(WL-VO0).GT.DWM) GOTO 130 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) c RZM=DSQRT(RRD(1)*RRD(1)+RRD(2)*RRD(2)+RRD(3)*RRD(3)) RZM=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) R1(1)=WL/PI180 R1(6)=ETT DWM=DABS(WL-VO0) C>>> 130 CONTINUE C>>> IF(DT.LT.DT05) GOTO 150 C>>> TA=R1(6)-DT C>>> TB=R1(6)+DT C>>> DT=0.1D-1*DT C>>> GOTO 120 C 150 QWV=FB*PI180 Ct WRITE(*,*) 'R1(6) = ',R1(6) Ct WLOUT=R1(1)+1.80D2 Ct IF(WLOUT.GE.3.60D2) WLOUT=WLOUT-3.60D2 Ct WRITE(*,*) 'JPL eph., lam = ',WLOUT Cx???? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) EB1(1)=RZM*(1.0D0+EA(2)*DCOS(QWV))/(1.0D0+EA(2)) IF(EB1(1).LT.0.0D0) GOTO 200 DD1=DABS(EA(1)-EB1(1)) R1(1)=R1(1)-0.18D3 IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3 VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS CALL RADIANT(EB1,FB,VXZ,VYR,VZR,R1) C 200 CONTINUE IF(EA(3).GT.0.18D3) GOTO 205 FB=0.36D3-EA(3) IF(FB.LT.1.800000001D2.AND.EA(2).GT.0.9999999D0) GOTO 300 VO0=EA(4) GOTO 210 205 FB=0.54D3-EA(3) VO0=EA(4)+0.18D3 IF(VO0.GE.0.36D3) VO0=VO0-0.36D3 210 VO0=VO0*PI180 DWM=2.0D0*PI2 TA=TMIN TB=TMAX DT=DT0 C 220 CONTINUE C>>> DO 230 ETT=TA,TB,DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 C>>> IF(DABS(WL-VO0).GT.DWM) GOTO 230 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) CX EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZM=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) R2(1)=WL/PI180 R2(6)=ETT DWM=DABS(WL-VO0) C>>> 230 CONTINUE C>>> IF(DT.LT.DT05) GOTO 250 C>>> TA=R2(6)-DT C>>> TB=R2(6)+DT C>>> DT=0.1D-1*DT C>>> GOTO 220 C 250 QWV=FB*PI180 Cx ???? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) EB2(1)=RZM*(1.0D0+EA(2)*DCOS(QWV))/(1.0D0+EA(2)) IF(EB2(1).LT.0.0D0) GOTO 300 DD2=DABS(EA(1)-EB2(1)) R2(1)=R2(1)-0.18D3 IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3 VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS CALL RADIANT(EB2,FB,VXZ,VYR,VZR,R2) C 300 CONTINUE RETURN C END C --------------------------------------------------------------------- SUBROUTINE BMETH(ETT,EA,R1,EB1,DD1,R2,EB2,DD2) C --------------------------------------------------------------------- C - (B) VARIATION OF PERIHELION DISTANCE AND EXCENTRICITY C (SVOREN J., NESLUSAN L., PORUBCAN V.: 1993, CONTRIB. ASTRON. OBS. C SKALNATE PLESO 23, 23.) C C "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE C POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 1950.0. NEXT C INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL- C DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2, C 3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS C THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D- C DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC C OF THE PARENT BODY ORBIT, RESPECTIVELY. C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) INTEGER I,IM,IYR,J DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7) Cx ,EZ(6) C PI=4.0D0*DATAN(1.0D0) PI2=2.0D0*PI PI180=PI/0.18D3 EPS=2.3439291D1 CEPS=DCOS(PI180*EPS) SEPS=DSIN(PI180*EPS) NTARG=3 NCENT=11 C DO 50 J=3,5 EB1(J)=EA(J) EB2(J)=EA(J) 50 CONTINUE DT0=(TMAX-TMIN)/3.6525D2 DT05=DT0/2.0D0 DD1=0.999D99 DD2=0.999D99 C IF(EA(3).GT.0.18D3) GOTO 100 FB=0.18D3-EA(3) IF(FB.GT.1.799999999D2.AND.EA(2).GT.0.9999999D0) GOTO 200 VO0=EA(4)+0.18D3 IF(VO0.GE.0.36D3) VO0=VO0-0.36D3 GOTO 110 100 FB=0.36D3-EA(3) VO0=EA(4) 110 VO0=VO0*PI180 DWM=2.0D0*PI2 TA=TMIN TB=TMAX DT=DT0 C 120 CONTINUE C>>> DO 130 ETT=TA,TB,DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 C>>> IF(DABS(WL-VO0).GT.DWM) GOTO 130 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZM=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R1(1)=WL/PI180 R1(6)=ETT DWM=DABS(WL-VO0) C>>> 130 CONTINUE C>>> IF(DT.LT.DT05) GOTO 150 C>>> TA=R1(6)-DT C>>> TB=R1(6)+DT C>>> DT=0.1D-1*DT C>>> GOTO 120 C 150 DCMIN=0.999D99 CF=DCOS(FB*PI180) Cx ??? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) C > DO 160 EC=0.0D0,5.0D0,0.1D0 DO 160 IEC=0,50 EC=IEC/1.0D1 DE=EC-EA(2) QX=RZM*(1.0D0+EC*CF)/(1.0D0+EC) DQ=EA(1)-QX DDX=DSQRT(DQ*DQ+DE*DE) IF(DDX.GE.DCMIN) GOTO 160 DCMIN=DDX EB1(2)=EC EB1(1)=QX 160 CONTINUE C ED=EB1(2)-0.2D0 IF(ED.LT.0.0D0) ED=0.0D0 EU=EB1(2)+0.2D0 C > DO 170 EC=ED,EU,0.1D-3 IEU=(EU+1.0D-11-ED)/0.1D-3 DO 170 IEC=0,IEU EC=ED+IEC*0.1D-3 DE=EC-EA(2) QX=RZM*(1.0D0+EC*CF)/(1.0D0+EC) DQ=EA(1)-QX DDX=DSQRT(DQ*DQ+DE*DE) IF(DDX.GE.DCMIN) GOTO 170 DCMIN=DDX EB1(2)=EC EB1(1)=QX 170 CONTINUE R1(1)=R1(1)-0.18D3 IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3 CALL RADIANT(EB1,FB,VXZ,VYR,VZR,R1) DD1=DCMIN C 200 CONTINUE IF(EA(3).GT.0.18D3) GOTO 205 FB=0.36D3-EA(3) IF(FB.LT.1.800000001D2.AND.EA(2).GT.0.9999999D0) GOTO 300 VO0=EA(4) GOTO 210 205 FB=0.54D3-EA(3) VO0=EA(4)+0.18D3 IF(VO0.GE.0.36D3) VO0=VO0-0.36D3 210 VO0=VO0*PI180 DWM=2.0D0*PI2 TA=TMIN TB=TMAX DT=DT0 C 220 CONTINUE C>>> DO 230 ETT=TA,TB,DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 C>>> IF(DABS(WL-VO0).GT.DWM) GOTO 230 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZM=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R2(1)=WL/PI180 R2(6)=ETT DWM=DABS(WL-VO0) C>>> 230 CONTINUE C>>> IF(DT.LT.DT05) GOTO 250 C>>> TA=R2(6)-DT C>>> TB=R2(6)+DT C>>> DT=0.1D-1*DT C>>> GOTO 220 C 250 DCMIN=0.999D99 CF=DCOS(FB*PI180) Cx ??? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) C > DO 260 EC=0.0D0,5.0D0,0.1D0 DO 260 IEC=0,50 EC=IEC/1.0D1 DE=EC-EA(2) QX=RZM*(1.0D0+EC*CF)/(1.0D0+EC) DQ=EA(1)-QX DDX=DSQRT(DQ*DQ+DE*DE) IF(DDX.GE.DCMIN) GOTO 260 DCMIN=DDX EB2(2)=EC EB2(1)=QX 260 CONTINUE C ED=EB2(2)-0.2D0 IF(ED.LT.0.0D0) ED=0.0D0 EU=EB2(2)+0.2D0 C > DO 270 EC=ED,EU,0.1D-3 IEU=(EU+1.0D-11-ED)/0.1D-3 DO 270 IEC=0,IEU EC=ED+IEC*0.1D-3 DE=EC-EA(2) QX=RZM*(1.0D0+EC*CF)/(1.0D0+EC) DQ=EA(1)-QX DDX=DSQRT(DQ*DQ+DE*DE) IF(DDX.GE.DCMIN) GOTO 270 DCMIN=DDX EB2(2)=EC EB2(1)=QX 270 CONTINUE R2(1)=R2(1)-0.18D3 IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3 CALL RADIANT(EB2,FB,VXZ,VYR,VZR,R2) DD2=DCMIN C 300 CONTINUE RETURN C END C --------------------------------------------------------------------- SUBROUTINE WMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2) C --------------------------------------------------------------------- C - (W) ROTATION OF THE LINE OF APSIDES C (STEEL D.I., BAGGALEY J.W.: 1985, MON. NOT. R. ASTRON. SOC. 212, C 817.) C C "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE C POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 1950.0. NEXT C INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL- C DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2, C 3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS C THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D- C DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC C OF THE PARENT BODY ORBIT, RESPECTIVELY. C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) INTEGER I,IM,IYR,J DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7) Cx ,EZ(6) C PI=4.0D0*DATAN(1.0D0) PI180=PI/0.18D3 PI05=PI/2.0D0 PI2=2.0D0*PI EPS=2.3439291D1 CEPS=DCOS(PI180*EPS) SEPS=DSIN(PI180*EPS) NTARG=3 NCENT=11 C DO 10 J=1,5 EB1(J)=EA(J) EB2(J)=EA(J) 10 CONTINUE IF(EA(2).GT.0.999D0) GOTO 40 IF(EA(1).GT.1.1D0) GOTO 20 IF(EA(1)*(1.0D0+EA(2))/(1.0D0-EA(2)).LT.0.9D0) GOTO 20 IF(EA(2).LT.1.0D-12) GOTO 20 GOTO 40 20 DD1=0.999D99 DD2=0.999D99 RETURN C 40 DT0=(TMAX-TMIN)/3.6525D2 DT05=DT0/2.0D0 C IF(EA(3).GT.0.18D3) GOTO 100 VO0=EA(4)+0.18D3 IF(VO0.GE.0.36D3) VO0=VO0-0.36D3 GOTO 110 100 VO0=EA(4) 110 VO0=VO0*PI180 DWM=2.0D0*PI2 TA=TMIN TB=TMAX DT=DT0 C 120 CONTINUE C > DO 130 ETT=TA,TB,DT IETTU=(TB+1.0D-11-TA)/DT DO 130 IETT=0,IETTU ETT=TA+IETT*DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 IF(DABS(WL-VO0).GT.DWM) GOTO 130 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZM=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R1(1)=WL/PI180 R1(6)=ETT DWM=DABS(WL-VO0) 130 CONTINUE IF(DT.LT.DT05) GOTO 150 TA=R1(6)-DT TB=R1(6)+DT DT=0.1D-1*DT GOTO 120 C 150 CONTINUE Cx ???? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) CFB=(EA(1)*(1.0D0+EA(2))/RZM-1.0D0)/EA(2) DD1=0.999D99 IF(DABS(CFB).GT.1.000001D0) GOTO 200 IF(CFB.LT.-1.0D0) CFB=-1.0D0 IF(CFB.GT.1.0D0) CFB=1.0D0 IF(DABS(CFB).LT.1.0D-12) GOTO 160 SFB=DSQRT(1.0D0-CFB*CFB) FB=DATAN(SFB/CFB) IF(FB.LT.0.0D0) FB=FB+PI GOTO 170 160 FB=PI05 170 FB=FB/PI180 EB1(3)=0.18D3-FB IF(EA(3).GT.0.18D3) EB1(3)=0.36D3-FB CALL DCRIT(EA,EB1,DD1) R1(1)=R1(1)-0.18D3 IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3 CALL RADIANT(EB1,FB,VXZ,VYR,VZR,R1) C 200 CONTINUE IF(EA(3).GT.0.18D3) GOTO 205 VO0=EA(4) GOTO 210 205 VO0=EA(4)+0.18D3 IF(VO0.GE.0.36D3) VO0=VO0-0.36D3 210 VO0=VO0*PI180 DWM=2.0D0*PI2 TA=TMIN TB=TMAX DT=DT0 C 220 CONTINUE C > DO 230 ETT=TA,TB,DT IETTU=(TB+1.0D-11-TA)/DT DO 230 IETT=0,IETTU ETT=TA+IETT*DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 IF(DABS(WL-VO0).GT.DWM) GOTO 230 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZM=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R2(1)=WL/PI180 R2(6)=ETT DWM=DABS(WL-VO0) 230 CONTINUE IF(DT.LT.DT05) GOTO 250 TA=R2(6)-DT TB=R2(6)+DT DT=0.1D-1*DT GOTO 220 C 250 CONTINUE Cx ???? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) CFB=(EA(1)*(1.0D0+EA(2))/RZM-1.0D0)/EA(2) DD2=0.999D99 IF(DABS(CFB).GT.1.000001D0) GOTO 300 IF(CFB.LT.-1.0D0) CFB=-1.0D0 IF(CFB.GT.1.0D0) CFB=1.0D0 IF(DABS(CFB).LT.1.0D-12) GOTO 260 SFB=DSQRT(1.0D0-CFB*CFB) FB=DATAN(SFB/CFB) IF(FB.LT.0.0D0) FB=FB+PI GOTO 270 260 FB=PI05 270 FB=0.36D3-FB/PI180 EB2(3)=0.36D3-FB IF(EA(3).GT.0.18D3) EB2(3)=0.54D3-FB CALL DCRIT(EA,EB2,DD2) R2(1)=R2(1)-0.18D3 IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3 CALL RADIANT(EB2,FB,VXZ,VYR,VZR,R2) C 300 CONTINUE RETURN C END C --------------------------------------------------------------------- SUBROUTINE AMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2) C --------------------------------------------------------------------- C - (A) ROTATION AROUND THE LINE OF APSIDES C (SVOREN J., NESLUSAN L., PORUBCAN V.: 1993, CONTRIB. ASTRON. OBS. C SKALNATE PLESO 23, 23.) C C "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE C POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 1950.0. NEXT C INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL- C DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2, C 3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS C THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D- C DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC C OF THE PARENT BODY ORBIT, RESPECTIVELY. C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) INTEGER I,IM,IYR,J DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7) Cx ,EZ(6) C PI=4.0D0*DATAN(1.0D0) PI05=PI/2.0D0 PI2=2.0D0*PI PI180=PI/0.18D3 EPS=2.3439291D1 CEPS=DCOS(PI180*EPS) SEPS=DSIN(PI180*EPS) NTARG=3 NCENT=11 C DD1=0.999D99 DD2=0.999D99 IF(EA(2).LT.1.0D-12) GOTO 5 IF(EA(2).GT.0.999D0) GOTO 10 IF(EA(1).LT.1.1D0) GOTO 10 IF(EA(1)*(1.0D0+EA(2))/(1.0D0-EA(2)).GT.0.9D0) GOTO 10 5 CONTINUE RETURN C 10 EB1(1)=EA(1) EB1(2)=EA(2) EB2(1)=EA(1) EB2(2)=EA(2) CSO=DCOS(PI180*EA(3)) SSO=DSIN(PI180*EA(3)) CHO=DCOS(PI180*EA(4)) SHO=DSIN(PI180*EA(4)) CSK=DCOS(PI180*EA(5)) SSK=DSIN(PI180*EA(5)) XQ=CSO*CHO-SSO*CSK*SHO YQ=CSO*SHO+SSO*CSK*CHO ZQ=SSO*SSK RQ=DSQRT(XQ*XQ+YQ*YQ+ZQ*ZQ) IF(DABS(RQ-1.0D0).GT.1.0D-12) WRITE(*,11) 11 FORMAT(/,' ERROR IN SUBROUTINE `AMETH` (RQ.NE.1).') IF(DABS(RQ-1.0D0).GT.1.0D-12) STOP C DT0=(TMAX-TMIN)/3.6525D2 DT05=DT0/2.0D0 C C > DO 13 ETT=TMIN,TMAX,DT05 IETTU=(TMAX+1.0D-11-TMIN)/DT05 DO 13 IETT=0,IETTUD ETT=TMIN+IETT*DT05 c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) Cx XZ=RRD(1) Cx YZ=RRD(2)*CEPS+RRD(3)*SEPS Cx ZZ=-RRD(2)*SEPS+RRD(3)*CEPS Cx VXZ=RRD(4) Cx VYZ=RRD(5)*CEPS+RRD(6)*SEPS Cx VZZ=-RRD(5)*SEPS+RRD(6)*CEPS Cx CALL ELEMS(ETT,XZ,YZ,VZ,VXZ,VYZ,VZZ,EZ) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) Cx RZ0=WEL(1)*(1.0D0-WEL(2)*WEL(2))/(1.0D0+WEL(2)*DCOS(WEL(7))) RZ0=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) CF0=(EA(1)*(1.0D0+EA(2))/RZ0-1.0D0)/EA(2) IF(DABS(CF0).LE.1.0D0) GOTO 15 13 CONTINUE GOTO 5 C 15 CONTINUE IF(DABS(CF0).LT.1.0D-9) GOTO 17 IF(DABS(CF0).GT.1.000000001D0) GOTO 200 IF(CF0.GT.1.0D0) CF0=1.0D0 IF(CF0.LT.-1.0D0) CF0=-1.0D0 SF0=DSQRT(1.0D0-CF0*CF0) F0=DATAN(SF0/CF0) IF(F0.LT.0.0D0) F0=F0+PI GOTO 20 17 F0=PI05 20 EB1(3)=PI-F0 IF(EA(3).GT.0.18D3) EB1(3)=PI2-F0 CSOB=DCOS(EB1(3)) SSOB=DSIN(EB1(3)) EB1(3)=EB1(3)/PI180 QKOB=SSOB*SSOB+CSOB*CSOB IF(DABS(QKOB-1.0D0).GT.1.0D-12) WRITE(*,30) 30 FORMAT(/,' ERROR IN SUBROUTINE `AMETH` (QKOB.NE.1).') IF(DABS(QKOB-1.0D0).GT.1.0D-12) STOP IF(DABS(SSOB).LT.1.0D-9.AND.DABS(SSO*SSK).GT.1.0D-9) GOTO 200 IF(DABS(SSOB).LT.1.0D-9) GOTO 40 SIB=SSO*SSK/SSOB IF(DABS(SIB).GT.1.0000000001D0) GOTO 200 GOTO 60 40 SIB=SSK 60 CONTINUE IF(SIB.GT.1.0D0) SIB=1.0D0 IF(SIB.LT.-1.0D0) SIB=-1.0D0 CIB=DSQRT(1.0D0-SIB*SIB) JIB=1 65 CONTINUE IF(DABS(CIB).LT.1.0D-12) EB1(5)=PI05 IF(DABS(CIB).LT.1.0D-12) GOTO 70 EB1(5)=DATAN(SIB/CIB) IF(CIB.LT.0.0D0) EB1(5)=PI+EB1(5) 70 CONTINUE EB1(5)=EB1(5)/PI180 WEVO=CSOB*CSOB+SSOB*SSOB*CIB*CIB C IF(DABS(CSOB).LT.1.0D-12.AND.DABS(CIB).LT.1.0D-12) DD1=0.999D99 IF(DABS(CSOB).LT.1.0D-12.AND.DABS(CIB).LT.1.0D-12) GOTO 200 QQS=(YQ*CSOB-XQ*SSOB*CIB)/WEVO QQC=(XQ*CSOB+YQ*SSOB*CIB)/WEVO QKONTR=QQS*QQS+QQC*QQC IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,80) 80 FORMAT(/,' ERROR IN SUBROUTINE `AMETH` ((SIN**2(EB(4)) + COS**(EB *(4))).GT.1).') IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP IF(DABS(QQC).LT.1.0D-12) EB1(4)=PI05 IF(DABS(QQC).LT.1.0D-12) GOTO 90 EB1(4)=DATAN(QQS/QQC) EB1(4)=DABS(EB1(4)) IF(QQC.GE.0.0D0.AND.QQS.LT.0.0D0) EB1(4)=PI2-EB1(4) IF(QQC.LT.0.0D0.AND.QQS.GE.0.0D0) EB1(4)=PI-EB1(4) IF(QQC.LT.0.0D0.AND.QQS.LT.0.0D0) EB1(4)=PI+EB1(4) 90 CONTINUE EB1(4)=EB1(4)/PI180 CALL DCRIT(EA,EB1,DD1) IF(JIB.EQ.2) GOTO 110 SOBO=EB1(3) HOBO=EB1(4) SKBO=EB1(5) DDO=DD1 100 JIB=2 CIB=-CIB GOTO 65 110 CONTINUE TA=TMIN TB=TMAX DT=DT0 IF(DD1.LT.DDO) GOTO 120 EB1(3)=SOBO EB1(4)=HOBO EB1(5)=SKBO DD1=DDO C 120 DWM=2.0D0*PI2 VO0=EB1(4)*PI180+PI IF(EB1(3).GT.0.18D3) VO0=EB1(4)*PI180 IF(VO0.GE.PI2) VO0=VO0-PI2 C > DO 130 ETT=TA,TB,DT IETTU=(TB+1.0D-11-TA)/DT DO 130 IETT=0,IETTU ETT=TA+IETT*DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 IF(DABS(WL-VO0).GT.DWM) GOTO 130 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZ=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R1(1)=WL/PI180 R1(6)=ETT DWM=DABS(WL-VO0) 130 CONTINUE IF(DT.LT.DT05) GOTO 150 TA=R1(6)-DT TB=R1(6)+DT DT=0.1D-1*DT GOTO 120 C 150 CONTINUE Cx ???? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) IF(DABS(RZ-RZ0).LT.1.0D-4) GOTO 160 RZ0=RZ CF0=(EA(1)*(1.0D0+EA(2))/RZ0-1.0D0)/EA(2) IF(DABS(CF0).LE.1.0D0) GOTO 15 GOTO 200 160 CONTINUE IF(DABS(CF0).LT.1.0D-9) GOTO 170 IF(DABS(CF0).GT.1.000000001D0) GOTO 200 IF(CF0.GT.1.0D0) CF0=1.0D0 IF(CF0.LT.-1.0D0) CF0=-1.0D0 SF0=DSQRT(1.0D0-CF0*CF0) F0=DATAN(SF0/CF0) IF(F0.LT.0.0D0) F0=F0+PI GOTO 180 170 F0=PI05 180 EB1(3)=(PI-F0)/PI180 IF(EA(3).GT.0.18D3) EB1(3)=(PI2-F0)/PI180 FB=F0/PI180 R1(1)=R1(1)-0.18D3 IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3 CALL RADIANT(EB1,FB,VXZ,VYR,VZR,R1) C C 2-ND ARC C 200 CONTINUE IF(DABS(CF0).LT.1.0D-9) GOTO 217 IF(DABS(CF0).GT.1.000000001D0) GOTO 400 IF(CF0.GT.1.0D0) CF0=1.0D0 IF(CF0.LT.-1.0D0) CF0=-1.0D0 SF0=DSQRT(1.0D0-CF0*CF0) F0=DATAN(SF0/CF0) IF(F0.LT.0.0D0) F0=F0+PI GOTO 220 217 F0=PI05 220 F0=PI2-F0 EB2(3)=PI2-F0 IF(EA(3).GT.0.18D3) EB2(3)=PI2+PI-F0 CSOB=DCOS(EB2(3)) SSOB=DSIN(EB2(3)) EB2(3)=EB2(3)/PI180 QKOB=SSOB*SSOB+CSOB*CSOB IF(DABS(QKOB-1.0D0).GT.1.0D-12) WRITE(*,230) 230 FORMAT(/,' ERROR IN SUBROUTINE `AMETH` (QKOB.NE.1).') IF(DABS(QKOB-1.0D0).GT.1.0D-12) STOP IF(DABS(SSOB).LT.1.0D-12.AND.DABS(SSO).GT.1.0D-9) GOTO 400 IF(DABS(SSOB).LT.1.0D-12) GOTO 245 SIB=SSO*SSK/SSOB IF(DABS(SIB).GT.1.0000000001D0) GOTO 400 GOTO 260 245 SIB=SSK 260 CONTINUE IF(SIB.GT.1.0D0) SIB=1.0D0 IF(SIB.LT.-1.0D0) SIB=-1.0D0 CIB=DSQRT(1.0D0-SIB*SIB) JIB=1 265 CONTINUE IF(DABS(CIB).LT.1.0D-12) EB2(5)=PI05 IF(DABS(CIB).LT.1.0D-12) GOTO 270 EB2(5)=DATAN(SIB/CIB) IF(CIB.LT.0.0D0) EB2(5)=PI+EB2(5) 270 CONTINUE EB2(5)=EB2(5)/PI180 WEVO=CSOB*CSOB+SSOB*SSOB*CIB*CIB C IF(DABS(CSOB).LT.1.0D-12.AND.DABS(CIB).LT.1.0D-12) DD2=0.999D99 IF(DABS(CSOB).LT.1.0D-12.AND.DABS(CIB).LT.1.0D-12) GOTO 400 QQS=(YQ*CSOB-XQ*SSOB*CIB)/WEVO QQC=(XQ*CSOB+YQ*SSOB*CIB)/WEVO QKONTR=QQS*QQS+QQC*QQC IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,280) 280 FORMAT(/,' ERROR IN SUBROUTINE `AMETH` ((SIN**2(EB(4)) + COS**(EB *(4))).GT.1).') IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP IF(DABS(QQC).LT.1.0D-12) EB2(4)=PI05 IF(DABS(QQC).LT.1.0D-12) GOTO 290 EB2(4)=DATAN(QQS/QQC) EB2(4)=DABS(EB2(4)) IF(QQC.GE.0.0D0.AND.QQS.LT.0.0D0) EB2(4)=PI2-EB2(4) IF(QQC.LT.0.0D0.AND.QQS.GE.0.0D0) EB2(4)=PI-EB2(4) IF(QQC.LT.0.0D0.AND.QQS.LT.0.0D0) EB2(4)=PI+EB2(4) 290 CONTINUE EB2(4)=EB2(4)/PI180 CALL DCRIT(EA,EB2,DD2) IF(JIB.EQ.2) GOTO 310 SOBO=EB2(3) HOBO=EB2(4) SKBO=EB2(5) DDO=DD2 300 JIB=2 CIB=-CIB GOTO 265 310 CONTINUE TA=TMIN TB=TMAX DT=DT0 IF(DD2.LT.DDO) GOTO 320 EB2(3)=SOBO EB2(4)=HOBO EB2(5)=SKBO DD2=DDO C 320 DWM=2.0D0*PI2 VO0=EB2(4)*PI180 IF(EB2(3).GT.0.18D3) VO0=EB2(4)*PI180+PI IF(VO0.GE.PI2) VO0=VO0-PI2 C > DO 330 ETT=TA,TB,DT IETTU=(TB+1.0D-11-TA)/DT DO 330 IETT=0,IETTU ETT=TA+IETT*DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 IF(DABS(WL-VO0).GT.DWM) GOTO 330 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZ=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R2(1)=WL/PI180 R2(6)=ETT DWM=DABS(WL-VO0) 330 CONTINUE IF(DT.LT.DT05) GOTO 350 TA=R2(6)-DT TB=R2(6)+DT DT=0.1D-1*DT GOTO 320 C 350 CONTINUE Cx ???? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) IF(DABS(RZ-RZ0).LT.1.0D-4) GOTO 360 RZ0=RZ CF0=(EA(1)*(1.0D0+EA(2))/RZ0-1.0D0)/EA(2) IF(DABS(CF0).LE.1.0D0) GOTO 200 GOTO 400 360 CONTINUE IF(DABS(CF0).LT.1.0D-9) GOTO 370 IF(DABS(CF0).GT.1.000000001D0) GOTO 400 IF(CF0.GT.1.0D0) CF0=1.0D0 IF(CF0.LT.-1.0D0) CF0=-1.0D0 SF0=DSQRT(1.0D0-CF0*CF0) F0=DATAN(SF0/CF0) IF(F0.LT.0.0D0) F0=F0+PI GOTO 380 370 F0=PI05 380 F0=PI2-F0 EB2(3)=(PI2-F0)/PI180 IF(EA(3).GT.0.18D3) EB2(3)=(PI2+PI-F0)/PI180 FB=F0/PI180 R2(1)=R2(1)-0.18D3 IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3 CALL RADIANT(EB2,FB,VXZ,VYR,VZR,R2) C 400 CONTINUE RETURN C END C --------------------------------------------------------------------- SUBROUTINE HMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2) C --------------------------------------------------------------------- C - (H) OMEGA-ADJUSTMENT C (HASEGAWA I.: 1990, PUBL. ASTRON. SOC. JAPAN 42, 175.) C C "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE C POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 1950.0. NEXT C INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL- C DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2, C 3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS C THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D- C DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC C OF THE PARENT BODY ORBIT, RESPECTIVELY. C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) INTEGER I,IM,IYR,J DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7) Cx ,EZ(6) C PI=4.0D0*DATAN(1.0D0) PI180=PI/0.18D3 PI05=PI/2.0D0 PI2=2.0D0*PI EPS=2.3439291D1 CEPS=DCOS(PI180*EPS) SEPS=DSIN(PI180*EPS) NTARG=3 NCENT=11 C GOTO 10 5 DD1=0.999D99 DD2=0.999D99 RETURN 10 EB1(1)=EA(1) EB1(2)=EA(2) EB2(1)=EA(1) EB2(2)=EA(2) CSK=DCOS(PI180*EA(5)) SSK=DSIN(PI180*EA(5)) C DT0=(TMAX-TMIN)/3.6525D2 DT05=DT0/2.0D0 C C > DO 13 ETT=TMIN,TMAX,DT05 IETTU=(TMAX+1.0D-11-TMIN)/DT05 DO 13 IETT=0,IETTU ETT=TMIN+IETT*DT05 c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) RZ0=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) CF0=(EA(1)*(1.0D0+EA(2))/RZ0-1.0D0)/EA(2) IF(DABS(CF0).LE.1.0D0) GOTO 15 13 CONTINUE GOTO 5 C 15 CONTINUE IF(DABS(CF0).LT.1.0D-9) GOTO 17 SF0=DSQRT(1.0D0-CF0*CF0) F0=DATAN(SF0/CF0) IF(F0.LT.0.0D0) F0=F0+PI GOTO 20 17 F0=PI05 20 U=EA(3)*PI180+F0 IF(U.GE.PI2) U=U-PI2 CU=DCOS(U) SU=DSIN(U) IF(DABS(CU).LT.1.0D-12) HLVO=PI05 IF(DABS(CU).LT.1.0D-12) GOTO 30 HLVO=DATAN(SU*CSK/CU) HLVO=DABS(HLVO) 30 CONTINUE IF(SU*CSK.GE.0.0D0.AND.CU.LT.0.0D0) HLVO=PI-HLVO IF(SU*CSK.LT.0.0D0.AND.CU.GE.0.0D0) HLVO=PI2-HLVO IF(SU*CSK.LT.0.0D0.AND.CU.LT.0.0D0) HLVO=PI+HLVO HL=HLVO+EA(4)*PI180 IF(HL.GE.PI2) HL=HL-PI2 C TA=TMIN TB=TMAX DT=DT0 DWM=2.0D0*PI2 C 40 CONTINUE C > DO 50 ETT=TA,TB,DT IETTU=(TB+1.0D-11-TA)/DT DO 50 IETT=0,IETTU ETT=TA+IETT*DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 IF(DABS(WL-HL).GT.DWM) GOTO 50 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZ=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R1(1)=WL/PI180 R1(6)=ETT DWM=DABS(WL-HL) 50 CONTINUE IF(DT.LT.DT05) GOTO 60 TA=R1(6)-DT TB=R1(6)+DT DT=0.1D-1*DT GOTO 40 C 60 CONTINUE Cx ???? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) CF0=(EA(1)*(1.0D0+EA(2))/RZ-1.0D0)/EA(2) IF(DABS(CF0).GT.1.000000001D0) GOTO 5 IF(DABS(RZ-RZ0).LT.1.0D-4) GOTO 65 RZ0=RZ GOTO 15 65 CONTINUE IF(DABS(CF0).LT.1.0D-9) GOTO 67 IF(CF0.GT.1.0D0) CF0=1.0D0 IF(CF0.LT.-1.0D0) CF0=-1.0D0 SF0=DSQRT(1.0D0-CF0*CF0) F0=DATAN(SF0/CF0) IF(F0.LT.0.0D0) F0=F0+PI GOTO 70 67 F0=PI05 C 70 CONTINUE EB1(3)=PI2-F0 EB1(4)=HL U=U/PI180 IF(U.LE.9.0D1.OR.U.GT.2.7D2) GOTO 80 EB1(3)=EB1(3)-PI IF(EB1(3).LT.0.0D0) EB1(3)=EB1(3)+PI2 EB1(4)=EB1(4)+PI IF(EB1(4).GE.PI2) EB1(4)=EB1(4)-PI2 80 EB1(3)=EB1(3)/PI180 EB1(4)=EB1(4)/PI180 C SIB=SSK*DCOS(HLVO) IF(DABS(SIB).GT.1.0000000001D0) GOTO 90 GOTO 100 90 CONTINUE WRITE(*,95) 95 FORMAT(/,' ERROR IN SUBROUTINE `HMETH` (ABS(SIN(EB(5))).GT.1).') STOP 100 CONTINUE IF(SIB.GT.1.0D0) SIB=1.0D0 IF(SIB.LT.-1.0D0) SIB=-1.0D0 CIB=DSQRT(1.0D0-SIB*SIB) IF(DABS(CIB).LT.1.0D-12) EB1(5)=PI05 IF(DABS(CIB).LT.1.0D-12) GOTO 110 EB1(5)=DATAN(SIB/CIB) IF(EB1(5).LT.0.0D0) EB1(5)=PI+EB1(5) 110 EB1(5)=EB1(5)/PI180 IF(EA(5).GT.9.0D1.AND.EB1(5).LT.9.0D1) EB1(5)=1.8D2-EB1(5) IF(EA(5).LT.9.0D1.AND.EB1(5).GT.9.0D1) EB1(5)=1.8D2-EB1(5) C CALL DCRIT(EA,EB1,DD1) FB=F0/PI180 R1(1)=R1(1)-0.18D3 IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3 CALL RADIANT(EB1,FB,VXZ,VYR,VZR,R1) C C 2-ND ARC C 200 CONTINUE IF(DABS(CF0).LT.1.0D-9) GOTO 217 SF0=DSQRT(1.0D0-CF0*CF0) F0=DATAN(SF0/CF0) IF(F0.LT.0.0D0) F0=F0+PI F0=PI2-F0 GOTO 220 217 F0=PI2-PI05 220 U=EA(3)*PI180+F0 IF(U.GE.PI2) U=U-PI2 CU=DCOS(U) SU=DSIN(U) IF(DABS(CU).LT.1.0D-12) HLVO=PI05 IF(DABS(CU).LT.1.0D-12) GOTO 230 HLVO=DATAN(SU*CSK/CU) HLVO=DABS(HLVO) 230 CONTINUE IF(SU*CSK.GE.0.0D0.AND.CU.LT.0.0D0) HLVO=PI-HLVO IF(SU*CSK.LT.0.0D0.AND.CU.GE.0.0D0) HLVO=PI2-HLVO IF(SU*CSK.LT.0.0D0.AND.CU.LT.0.0D0) HLVO=PI+HLVO HL=HLVO+EA(4)*PI180 IF(HL.GE.PI2) HL=HL-PI2 C TA=TMIN TB=TMAX DT=DT0 DWM=2.0D0*PI2 C 240 CONTINUE C > DO 250 ETT=TA,TB,DT IETTU=(TB+1.0D-11-TA)/DT DO 250 IETT=0,IETTU ETT=TA+IETT*DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 IF(DABS(WL-HL).GT.DWM) GOTO 250 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZ=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R2(1)=WL/PI180 R2(6)=ETT DWM=DABS(WL-HL) 250 CONTINUE IF(DT.LT.DT05) GOTO 260 TA=R2(6)-DT TB=R2(6)+DT DT=0.1D-1*DT GOTO 240 C 260 CONTINUE Cx ????? IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM)) CF0=(EA(1)*(1.0D0+EA(2))/RZ-1.0D0)/EA(2) IF(DABS(CF0).GT.1.000000001D0) GOTO 5 IF(DABS(RZ-RZ0).LT.1.0D-4) GOTO 265 RZ0=RZ GOTO 200 265 CONTINUE IF(DABS(CF0).LT.1.0D-9) GOTO 267 IF(CF0.GT.1.0D0) CF0=1.0D0 IF(CF0.LT.-1.0D0) CF0=-1.0D0 SF0=DSQRT(1.0D0-CF0*CF0) F0=DATAN(SF0/CF0) IF(F0.LT.0.0D0) F0=F0+PI F0=PI2-F0 GOTO 270 267 F0=PI2-PI05 C 270 CONTINUE EB2(3)=PI2-F0 EB2(4)=HL U=U/PI180 IF(U.LE.9.0D1.OR.U.GT.2.7D2) GOTO 280 EB2(3)=EB2(3)-PI IF(EB2(3).LT.0.0D0) EB2(3)=EB2(3)+PI2 EB2(4)=EB2(4)+PI IF(EB2(4).GE.PI2) EB2(4)=EB2(4)-PI2 280 EB2(3)=EB2(3)/PI180 EB2(4)=EB2(4)/PI180 C SIB=SSK*DCOS(HLVO) IF(DABS(SIB).GT.1.0000000001D0) GOTO 290 GOTO 300 290 CONTINUE WRITE(*,95) STOP 300 CONTINUE IF(SIB.GT.1.0D0) SIB=1.0D0 IF(SIB.LT.-1.0D0) SIB=-1.0D0 CIB=DSQRT(1.0D0-SIB*SIB) IF(DABS(CIB).LT.1.0D-12) EB2(5)=PI05 IF(DABS(CIB).LT.1.0D-12) GOTO 310 EB2(5)=DATAN(SIB/CIB) IF(EB2(5).LT.0.0D0) EB2(5)=PI+EB2(5) 310 EB2(5)=EB2(5)/PI180 IF(EA(5).GT.9.0D1.AND.EB2(5).LT.9.0D1) EB2(5)=1.8D2-EB2(5) IF(EA(5).LT.9.0D1.AND.EB2(5).GT.9.0D1) EB2(5)=1.8D2-EB2(5) C CALL DCRIT(EA,EB2,DD2) FB=F0/PI180 R2(1)=R2(1)-0.18D3 IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3 CALL RADIANT(EB2,FB,VXZ,VYR,VZR,R2) C 400 CONTINUE RETURN C END C --------------------------------------------------------------------- SUBROUTINE PMETH(WM1,WM2,EA,R1,EB1,DD1,R2,EB2,DD2) C --------------------------------------------------------------------- C - (P) PORTER'S METHOD C (PORTER J.G.: 1952, COMETS AND METEOR STREAMS, CHAPMAN AND HALL C LTD., LONDON.) C C INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL- C DS OF VALUES "R1(I)", "R2(I)", "WM1(I)", WM2(I)", "EB1(I)", AND C "EB2(I)", I = 1, 2, 3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE C MAIN PROGRAM. AS THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" C AND "DD2" OF D-DISCRIMINANT CORRESPONDING TO THE POST- AND PRE- C PERIHELION ARC OF THE PARENT BODY ORBIT, RESPECTIVELY. C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) INTEGER I,IM,IYR,J DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WM1(6),WM2(6),WEL(7) Cx DIMENSION EZ(6) C PI=4.0D0*DATAN(1.0D0) PI05=PI/2.0D0 PI2=2.0D0*PI PI180=PI/0.18D3 GAUSS=.17202098950D-1 GS2=GAUSS*GAUSS EPS=2.3439291D1 CEPS=DCOS(PI180*EPS) SEPS=DSIN(PI180*EPS) NTARG=3 NCENT=11 C CHO=DCOS(EA(4)*PI180) SHO=DSIN(EA(4)*PI180) CSK=DCOS(EA(5)*PI180) SSK=DSIN(EA(5)*PI180) C DD1=0.999D99 IF(WM1(2).GT.1.0D2.OR.WM1(5).LT.0.0D0) GOTO 200 CVK=DCOS(WM1(5)*PI180) SVK=DSIN(WM1(5)*PI180) WVFI=DSQRT(1.0D0+EA(2)*EA(2)+2.0D0*EA(2)*CVK) CHFI=-EA(2)*SVK/WVFI SHFI=(1.0D0+EA(2)*CVK)/WVFI QKONTR=SHFI*SHFI+CHFI*CHFI IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,5) 5 FORMAT(/,' ERROR IN SUBROUTINE `PMETH` ((SIN**2(HFI) + COS**2(VFI *).NE.1).') IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP IF(SHFI.LT.-1.0D-12) WRITE(*,7) 7 FORMAT(/,' ERROR IN SUBROUTINE `PMETH` (ANGLE `HFI` > PI).') IF(SHFI.LT.-1.0D-12) STOP IF(DABS(CHFI).LT.1.0D-12) HFI=PI05 IF(DABS(CHFI).LT.1.0D-12) GOTO 10 HFI=DATAN(SHFI/CHFI) IF(HFI.LT.0.0D0) HFI=PI+HFI 10 TMN1=WM1(1) c CALL PLEPH(TMN1,NTARG,NCENT,RRD) T20=(TMN1-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) RZ=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R1(1)=WL/PI180 R1(6)=TMN1 Cx RZ=DSQRT(X*X+Y*Y+Z*Z) IF(EA(1).GT.RZ) GOTO 20 IF(EA(2).GE.0.9999999D0) GOTO 30 QQK=EA(1)*(1.0D0+EA(2))/(1.0D0-EA(2)) IF(QQK.GE.RZ) GOTO 30 AV=GAUSS*DSQRT(2.0D0/QQK-(1.0D0-EA(2))/EA(1)) GOTO 40 20 CONTINUE AV=GAUSS*DSQRT((1.0D0+EA(2))/EA(1)) GOTO 40 30 CONTINUE AV=GAUSS*DSQRT(2.0D0/RZ-(1.0D0-EA(2))/EA(1)) 40 AR3S=(WM1(5)+EA(3))*PI180-HFI C3S=DCOS(AR3S) S3S=DSIN(AR3S) VX=-AV*(C3S*CHO-S3S*CSK*SHO) VY=-AV*(C3S*SHO+S3S*CSK*CHO) VZ=-AV*S3S*SSK C WMEN=(XZ*VY-YZ*VX)*(XZ*VY-YZ*VX)+(YZ*VZ-ZZ*VY)*(YZ*VZ-ZZ*VY) WMEN=DSQRT(WMEN+(ZZ*VX-XZ*VZ)*(ZZ*VX-XZ*VZ)) CIB=(XZ*VY-YZ*VX)/WMEN SIB=DSQRT(1.0D0-CIB*CIB) IF(DABS(CIB).LT.1.0D-12) GOTO 50 EB1(5)=DATAN(SIB/CIB) IF(EB1(5).LT.0.0D0) EB1(5)=PI+EB1(5) GOTO 60 50 EB1(5)=PI05 60 WMEN=(YZ*VZ-ZZ*VY)*(YZ*VZ-ZZ*VY)+(ZZ*VX-XZ*VZ)*(ZZ*VX-XZ*VZ) WMEN=DSQRT(WMEN) CVO=(XZ*VZ-ZZ*VX)/WMEN SVO=(YZ*VZ-ZZ*VY)/WMEN IF(DABS(CVO).LT.1.0D-12) GOTO 70 EB1(4)=DATAN(SVO/CVO) EB1(4)=DABS(EB1(4)) GOTO 80 70 EB1(4)=PI05 80 CONTINUE IF(CVO.GE.0.0D0.AND.SVO.LT.0.0D0) EB1(4)=PI2-EB1(4) IF(CVO.LT.0.0D0.AND.SVO.GE.0.0D0) EB1(4)=PI-EB1(4) IF(CVO.LT.0.0D0.AND.SVO.LT.0.0D0) EB1(4)=PI+EB1(4) C AKI=2.0D0/RZ-AV*AV/GS2 EB1(2)=DSQRT(1.0D0-RZ*AKI*(2.0D0-RZ*AKI)*SHFI*SHFI) IF(DABS(EB1(2)-1.0D0).LT.1.0D-9) GOTO 90 EB1(1)=(1.0D0-EB1(2))/AKI GOTO 100 90 EB1(1)=RZ*SHFI*SHFI 100 CF=(EB1(1)*(1.0D0+EB1(2))/RZ-1.0D0)/EB1(2) IF(CF.GE.0.999999999D0) GOTO 115 IF(CF.LE.-0.999999999D0) GOTO 117 SF=DSQRT(1.0D0-CF*CF) IF(DABS(CF).LT.1.0D-12) GOTO 110 F=DATAN(SF/CF) IF(F.LT.0.0D0) F=PI+F GOTO 120 110 F=PI05 GOTO 120 115 F=0.0D0 GOTO 120 117 F=PI 120 CONTINUE IF(HFI.LT.PI05) F=PI2-F C C2S=(XZ*CVO+YZ*SVO)/RZ IF(DABS(CIB).LT.1.0D-9) GOTO 130 S2S=(YZ*CVO-XZ*SVO)/RZ/CIB GOTO 140 130 S2S=ZZ/RZ/SIB 140 CONTINUE IF(DABS(C2S).LT.1.0D-12) GOTO 150 U2S=DATAN(S2S/C2S) U2S=DABS(U2S) GOTO 160 150 U2S=PI05 160 CONTINUE IF(C2S.GE.0.0D0.AND.S2S.LT.0.0D0) U2S=PI2-U2S IF(C2S.LT.0.0D0.AND.S2S.GE.0.0D0) U2S=PI-U2S IF(C2S.LT.0.0D0.AND.S2S.LT.0.0D0) U2S=PI+U2S EB1(3)=U2S-F IF(EB1(3).LT.0.0D0) EB1(3)=EB1(3)+PI2 C EB1(3)=EB1(3)/PI180 EB1(4)=EB1(4)/PI180 EB1(5)=EB1(5)/PI180 CALL DCRIT(EA,EB1,DD1) FB=F/PI180 R1(1)=R1(1)-0.18D3 IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3 CALL RADIANT(EB1,FB,VXZ,VYR,VZR,R1) C C 2-ND ARC C 200 CONTINUE DD2=0.999D99 IF(WM2(2).GT.1.0D2.OR.WM2(5).GT.0.36D3) RETURN CVK=DCOS(WM2(5)*PI180) SVK=DSIN(WM2(5)*PI180) WVFI=DSQRT(1.0D0+EA(2)*EA(2)+2.0D0*EA(2)*CVK) CHFI=-EA(2)*SVK/WVFI SHFI=(1.0D0+EA(2)*CVK)/WVFI QKONTR=SHFI*SHFI+CHFI*CHFI IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,5) IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP IF(SHFI.LT.-1.0D-12) WRITE(*,7) IF(SHFI.LT.-1.0D-12) STOP IF(DABS(CHFI).LT.1.0D-12) HFI=PI05 IF(DABS(CHFI).LT.1.0D-12) GOTO 210 HFI=DATAN(SHFI/CHFI) IF(HFI.LT.0.0D0) HFI=PI+HFI 210 TMN2=WM2(1) c CALL PLEPH(TMN2,NTARG,NCENT,RRD) T20=(TMN2-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) IF(DABS(XZ).LT.1.0D-12) THEN WL=PI/2.0D0 IF(YZ.LT.0.0D0) WL=WL+PI ELSE WL=DATAN(YZ/XZ) WL=DABS(WL) IF(XZ.LT.0.0D0.AND.YZ.GE.0.0D0) WL=PI-WL IF(XZ.LT.0.0D0.AND.YZ.LT.0.0D0) WL=PI+WL IF(XZ.GE.0.0D0.AND.YZ.LT.0.0D0) WL=PI2-WL END IF C WL=WEL(7)+WEL(4) C IF(WL.GE.PI2) WL=WL-PI2 Cx EZ(1)=WEL(1) Cx EZ(2)=WEL(2) Cx EZ(3)=WEL(4)-WEL(5) Cx IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2 Cx EZ(4)=WEL(5) Cx EZ(5)=WEL(6) Cx FZM=WEL(7) Cx EZ(1)=EZ(1)/(1.0D0-EZ(2)) Cx FZM=EZ(6) VYR=VYZ*CEPS-VZZ*SEPS VZR=VYZ*SEPS+VZZ*CEPS R2(1)=WL/PI180 R2(6)=TMN2 RZ=DSQRT(XZ*XZ+YZ*YZ+ZZ*ZZ) IF(EA(1).GT.RZ) GOTO 220 IF(EA(2).GE.0.9999999D0) GOTO 230 QQK=EA(1)*(1.0D0+EA(2))/(1.0D0-EA(2)) IF(QQK.GE.RZ) GOTO 230 AV=GAUSS*DSQRT(2.0D0/QQK-(1.0D0-EA(2))/EA(1)) GOTO 240 220 CONTINUE AV=GAUSS*DSQRT((1.0D0+EA(2))/EA(1)) GOTO 240 230 CONTINUE AV=GAUSS*DSQRT(2.0D0/RZ-(1.0D0-EA(2))/EA(1)) 240 AR3S=(WM2(5)+EA(3))*PI180-HFI C3S=DCOS(AR3S) S3S=DSIN(AR3S) VX=-AV*(C3S*CHO-S3S*CSK*SHO) VY=-AV*(C3S*SHO+S3S*CSK*CHO) VZ=-AV*S3S*SSK C WMEN=(XZ*VY-YZ*VX)*(XZ*VY-YZ*VX)+(YZ*VZ-ZZ*VY)*(YZ*VZ-ZZ*VY) WMEN=DSQRT(WMEN+(ZZ*VX-XZ*VZ)*(ZZ*VX-XZ*VZ)) CIB=(XZ*VY-YZ*VX)/WMEN SIB=DSQRT(1.0D0-CIB*CIB) IF(DABS(CIB).LT.1.0D-12) GOTO 250 EB2(5)=DATAN(SIB/CIB) IF(EB2(5).LT.0.0D0) EB2(5)=PI+EB2(5) GOTO 260 250 EB2(5)=PI05 260 WMEN=(YZ*VZ-ZZ*VY)*(YZ*VZ-ZZ*VY)+(ZZ*VX-XZ*VZ)*(ZZ*VX-XZ*VZ) WMEN=DSQRT(WMEN) CVO=(XZ*VZ-ZZ*VX)/WMEN SVO=(YZ*VZ-ZZ*VY)/WMEN IF(DABS(CVO).LT.1.0D-12) GOTO 270 EB2(4)=DATAN(SVO/CVO) EB2(4)=DABS(EB2(4)) GOTO 280 270 EB2(4)=PI05 280 CONTINUE IF(CVO.GE.0.0D0.AND.SVO.LT.0.0D0) EB2(4)=PI2-EB2(4) IF(CVO.LT.0.0D0.AND.SVO.GE.0.0D0) EB2(4)=PI-EB2(4) IF(CVO.LT.0.0D0.AND.SVO.LT.0.0D0) EB2(4)=PI+EB2(4) C AKI=2.0D0/RZ-AV*AV/GS2 EB2(2)=DSQRT(1.0D0-RZ*AKI*(2.0D0-RZ*AKI)*SHFI*SHFI) IF(DABS(EB2(2)-1.0D0).LT.1.0D-9) GOTO 290 EB2(1)=(1.0D0-EB2(2))/AKI GOTO 300 290 EB2(1)=RZ*SHFI*SHFI 300 CF=(EB2(1)*(1.0D0+EB2(2))/RZ-1.0D0)/EB2(2) IF(CF.GE.0.999999999D0) GOTO 315 IF(CF.LE.-0.999999999D0) GOTO 317 SF=DSQRT(1.0D0-CF*CF) IF(DABS(CF).LT.1.0D-12) GOTO 310 F=DATAN(SF/CF) IF(F.LT.0.0D0) F=PI+F GOTO 320 310 F=PI05 GOTO 320 315 F=0.0D0 GOTO 320 317 F=PI 320 CONTINUE IF(HFI.LT.PI05) F=PI2-F C C2S=(XZ*CVO+YZ*SVO)/RZ IF(DABS(CIB).LT.1.0D-9) GOTO 330 S2S=(YZ*CVO-XZ*SVO)/RZ/CIB GOTO 340 330 S2S=ZZ/RZ/SIB 340 CONTINUE IF(DABS(C2S).LT.1.0D-12) GOTO 350 U2S=DATAN(S2S/C2S) U2S=DABS(U2S) GOTO 360 350 U2S=PI05 360 CONTINUE IF(C2S.GE.0.0D0.AND.S2S.LT.0.0D0) U2S=PI2-U2S IF(C2S.LT.0.0D0.AND.S2S.GE.0.0D0) U2S=PI-U2S IF(C2S.LT.0.0D0.AND.S2S.LT.0.0D0) U2S=PI+U2S EB2(3)=U2S-F IF(EB2(3).LT.0.0D0) EB2(3)=EB2(3)+PI2 C EB2(3)=EB2(3)/PI180 EB2(4)=EB2(4)/PI180 EB2(5)=EB2(5)/PI180 CALL DCRIT(EA,EB2,DD2) FB=F/PI180 R2(1)=R2(1)-0.18D3 IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3 CALL RADIANT(EB2,FB,VXZ,VYR,VZR,R2) C RETURN C END C --------------------------------------------------------------------- SUBROUTINE RADIANT(EB,FB,VXZ,VYZ,VZZ,RR) C --------------------------------------------------------------------- C - CALCULATION OF METEOR SHOWER RADIANT, IF THE FOLLOWING INPUT C QUANTITIES OF THE ORBIT CROSSING THE EARTH'S ORBIT ARE WELL-KNOWN: C C EB(1) - PERIHELION DISTANCE [AU] C EB(2) - ECCENTRICITY C EB(3) - ARGUMENT OF PERIHELION [DEGREES] C EB(4) - LONGITUDE OF NODE [DEGREES] C EB(5) - INCLINATION TO ECLIPTIC [DEGREES] C FB - TRUE ANOMALY OF THE POINT OF INTERSECTION OF CALCULATED C SHOWER AND EARTH ORBITS [DEGREES] C VXZ, VYZ, VZZ - HELIOCENTIC EQUATORIAL RECTANULAR COMPONENTS C IF THE EARTH [AU/DAY] C C OUTPUT: C THE QUANTITIES DESCRIBING THE RADIANT OF POTENTIAL METEOR SHOWER C RR(2) - RIGHT ASCENSION OF THE PREDICTED RADIANT [DEGREES] C RR(3) - DECLINATION OF THE RADIANT [DEGREES] C RR(4) - PREDICTED GEOCENTRIC VELOCITY OF METEORS [KM/S] C RR(5) - PREDICTED HELIOCENTRIC VELOCITY OF METEORS [KM/S] C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) DIMENSION RR(6),EB(5) Cx ,EZ(6) C PI=4.0D0*DATAN(1.0D0) PI05=PI/2.0D0 PI2=2.0D0*PI PI180=PI/0.18D3 GAUSS=0.17202098950D-1 AUKM=0.149597870D9 D24=8.6400D4 C EPS=2.344579D1......1950.0 EPS=2.3439291D1 C CEPS=DCOS(PI180*EPS) SEPS=DSIN(PI180*EPS) CHO=DCOS(PI180*EB(4)) SHO=DSIN(PI180*EB(4)) CSK=DCOS(PI180*EB(5)) SSK=DSIN(PI180*EB(5)) RK=EB(1)*(1.0D0+EB(2))/(1.0D0+EB(2)*DCOS(FB*PI180)) Cx CHZ=DCOS(EZ(4)) Cx SHZ=DSIN(EZ(4)) Cx CSZ=DCOS(EZ(5)) Cx SSZ=DSIN(EZ(5)) Cx RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZ)) C VHR=PI180*FB CVH=DCOS(VHR) SVH=DSIN(VHR) WHFI=DSQRT(1.0D0+EB(2)*EB(2)+2.0D0*EB(2)*CVH) CHFI=-EB(2)*SVH/WHFI SHFI=(1.0D0+EB(2)*CVH)/WHFI QKONTR=SHFI*SHFI+CHFI*CHFI IF(DABS(QKONTR-1.0D0).GT.1.0D-9) THEN WRITE(*,100) 100 FORMAT(/,' ERROR IN SUBROUTINE `RADIANT` (ABS(SIN**2(HFI) + COS** *2(HFI)).GT.1)') WRITE(*,*) 'QKONTR =',QKONTR ENDIF IF(DABS(QKONTR-1.0D0).GT.1.0D-9) STOP IF(DABS(CHFI).LT.1.0D-9) HFI=PI05 IF(DABS(CHFI).LT.1.0D-9) GOTO 120 HFI=DATAN(SHFI/CHFI) HFI=DABS(HFI) IF(CHFI.GE.0.0D0.AND.SHFI.LT.0.0D0) HFI=PI2-HFI IF(CHFI.LT.0.0D0.AND.SHFI.GE.0.0D0) HFI=PI-HFI IF(CHFI.LT.0.0D0.AND.SHFI.LT.0.0D0) HFI=PI+HFI 120 CONTINUE CFSO=DCOS(PI180*(FB+EB(3))) SFSO=DSIN(PI180*(FB+EB(3))) XK=RK*(CFSO*CHO-SFSO*CSK*SHO) YK=RK*(CFSO*SHO+SFSO*CSK*CHO) ZK=RK*SFSO*SSK C3=DCOS(PI180*(FB+EB(3))-HFI) S3=DSIN(PI180*(FB+EB(3))-HFI) AV=GAUSS*DSQRT(2.0D0/RK-(1.0D0-EB(2))/EB(1)) VXE=-AV*(C3*CHO-S3*CSK*SHO) VYE=-AV*(C3*SHO+S3*CSK*CHO) VZE=-AV*S3*SSK VE=DSQRT(VXE*VXE+VYE*VYE+VZE*VZE) IF(DABS(AV-VE).GT.1.0D-9) WRITE(*,140) 140 FORMAT(/,' ERROR IN SUBROUTINE `RADIANT` (AV.NE.VE).') IF(DABS(AV-VE).GT.1.0D-9) STOP VXK=VXE VYK=VYE*CEPS-VZE*SEPS VZK=VYE*SEPS+VZE*CEPS VH=DSQRT(VXK*VXK+VYK*VYK+VZK*VZK) IF(DABS(AV-VH).GT.1.0D-9) WRITE(*,160) 160 FORMAT(/,' ERROR IN SUBROUTINE `RADIANT` (AV.NE.VH).') IF(DABS(AV-VH).GT.1.0D-9) STOP C Cx CVH=DCOS(FZ) Cx SVH=DSIN(FZ) Cx WHFI=DSQRT(1.0D0+EZ(2)*EZ(2)+2.0D0*EZ(2)*CVH) Cx CHFI=-EZ(2)*SVH/WHFI Cx SHFI=(1.0D0+EZ(2)*CVH)/WHFI Cx QKONTR=SHFI*SHFI+CHFI*CHFI Cx IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,200) Cx 200 FORMAT(/,' ERROR IN SUBROUTINE `RADIANT` (ABS(SIN**2(HFI) + COS** Cx *2(HFI)).GT.1)') Cx IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP Cx IF(DABS(CHFI).LT.1.0D-12) HFI=PI05 Cx IF(DABS(CHFI).LT.1.0D-12) GOTO 220 Cx HFI=DATAN(SHFI/CHFI) Cx HFI=DABS(HFI) Cx IF(CHFI.GE.0.0D0.AND.SHFI.LT.0.0D0) HFI=PI2-HFI Cx IF(CHFI.LT.0.0D0.AND.SHFI.GE.0.0D0) HFI=PI-HFI Cx IF(CHFI.LT.0.0D0.AND.SHFI.LT.0.0D0) HFI=PI+HFI Cx 220 CONTINUE Cx CFSO=DCOS(FZ+EZ(3)) Cx SFSO=DSIN(FZ+EZ(3)) Cx XZ=RZ*(CFSO*CHZ-SFSO*CSZ*SHZ) Cx YZ=RZ*(CFSO*SHZ+SFSO*CSZ*CHZ) Cx ZZ=RZ*SFSO*SSZ Cx C3=DCOS(FZ+EZ(3)-HFI) Cx S3=DSIN(FZ+EZ(3)-HFI) Cx AVZ=GAUSS*DSQRT(2.0D0/RZ-1.0D0/EZ(1)) Cx VXE=-AVZ*(C3*CHZ-S3*CSZ*SHZ) Cx VYE=-AVZ*(C3*SHZ+S3*CSZ*CHZ) Cx VZE=-AVZ*S3*SSZ Cx VEZ=DSQRT(VXE*VXE+VYE*VYE+VZE*VZE) Cx IF(DABS(AVZ-VEZ).GT.1.0D-12) WRITE(*,240) Cx 240 FORMAT(/,' ERROR IN SUBROUTINE `RADIANT` (AVZ.NE.VEZ).') Cx IF(DABS(AVZ-VEZ).GT.1.0D-12) STOP Cx VXZ=VXE Cx VYZ=VYE*CEPS-VZE*SEPS Cx VZZ=VYE*SEPS+VZE*CEPS Cx VHZ=DSQRT(VXZ*VXZ+VYZ*VYZ+VZZ*VZZ) Cx IF(DABS(AVZ-VHZ).GT.1.0D-12) WRITE(*,260) Cx 260 FORMAT(/,' ERROR IN SUBROUTINE `RADIANT` (AVZ.NE.VHZ).') Cx IF(DABS(AVZ-VHZ).GT.1.0D-12) STOP C VX=-(VXK-VXZ) VY=-(VYK-VYZ) VZ=-(VZK-VZZ) VG=DSQRT(VX*VX+VY*VY+VZ*VZ) IF(DABS(VZ).LT.1.0D-9) D=0.0D0 IF(DABS(VZ).LT.1.0D-9) GOTO 310 HELVEC=DSQRT(1.0D0-VZ*VZ/VG/VG) D=DATAN(VZ/HELVEC/VG) 310 CONTINUE IF(DABS(D-PI05).LT.1.0D-9.OR.DABS(D+PI05).LT.1.0D-9) A=0.0D0 IF(DABS(D-PI05).LT.1.0D-9.OR.DABS(D+PI05).LT.1.0D-9) GOTO 340 CDVEC=DCOS(D) QQS=VY/VG/CDVEC QQC=VX/VG/CDVEC QKONTR=QQS*QQS+QQC*QQC IF(DABS(QKONTR-1.0D0).GT.1.0D-9) WRITE(*,320) 320 FORMAT(/,' ERROR IN SUBROUTINE `RADIANT` (SIN**2(A) + COS**2(A).G *T.1).') IF(DABS(QKONTR-1.0D0).GT.1.0D-9) STOP IF(DABS(QQC).LT.1.0D-9) A=PI05 IF(DABS(QQC).LT.1.0D-9) GOTO 340 A=DATAN(QQS/QQC) A=DABS(A) IF(QQC.GE.0.0D0.AND.QQS.LT.0.0D0) A=PI2-A IF(QQC.LT.0.0D0.AND.QQS.GE.0.0D0) A=PI-A IF(QQC.LT.0.0D0.AND.QQS.LT.0.0D0) A=PI+A 340 CONTINUE RR(2)=A/PI180 RR(3)=D/PI180 RR(4)=VG*AUKM/D24 RR(5)=AV*AUKM/D24 RETURN C END C --------------------------------------------------------------------- SUBROUTINE DCRIT(EA,EB,DD) C --------------------------------------------------------------------- C - COMPUTATION OF D-DISCRIMINANT OF THE SOUTHWORTH-HAWKINS' D-CRI- C TERION EVALUATING SIMILARITY OF TWO ORBITS (SOUTHWORTH R.B., C HAWKINS G.S.: 1963, SMITHSON. CONTR. ASTROPHYS. 7, 261.) C C INPUT: C EA(1), EB(1) - PERIHELION DISTANCES OF THE ORBITS [AU] C EA(2), EB(2) - ECCENTRICITIES C EA(3), EB(3) - ARGUMENTS OF PERIHELIA [DEGREES] C EA(4), EB(4) - LONGITUDES OF NODES [DEGREES] C EA(5), EB(5) - INCLINATIONS TO ECLIPTIC [DEGREES] C C OUTPUT: C DD - RESULTANT VALUE OF D-DISCRIMINANT C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8(A-H,O-Z) DIMENSION EA(6),EB(5) C PI=4.0D0*DATAN(1.0D0) PI180=PI/0.18D3 C DQ=EA(1)-EB(1) DE=EA(2)-EB(2) DSO=(EA(3)-EB(3))*PI180 DHO=(EA(4)-EB(4))*PI180 DSK=(EA(5)-EB(5))*PI180 SSKA=DSIN(EA(5)*PI180) SSKB=DSIN(EB(5)*PI180) POMSK=DSIN(DSK/2.0D0) POMHO=DSIN(DHO/2.0D0) DAI=4.0D0*POMSK*POMSK+SSKA*SSKB*4.0D0*POMHO*POMHO IF(DAI.GE.3.9999999D0.AND.DAI.LT.4.0000001D0) DAI=3.9999999D0 IF(DAI.LT.4.0D0) GOTO 700 WRITE(*,600) 600 FORMAT(/,' ERROR IN SUBROUTINE `DCRIT` (QUANTITY DAI.GE.4).') WRITE(*,*) DAI STOP 700 SECSK=1.0D0/SQRT(1.0D0-DAI/4.0D0) ARGS=DCOS((EA(5)+EB(5))*PI180/2.0D0)*POMHO*SECSK JSGAR=-1 IF(ARGS.GE.0.0D0) JSGAR=1 POM=DABS(DHO) JSIGN=2 IF(POM.GT.PI) JSIGN=-2 POM=DABS(ARGS) IF(POM.LT.1.0000001D0) GOTO 720 WRITE(*,620) 620 FORMAT(/,' ERROR IN SUBROUTINE `DCRIT` (ABS(ARGS).GT.1).') STOP 720 CONTINUE IF(POM.LT.1) GOTO 740 DBP=DSO+JSIGN*JSGAR*PI/2.0D0 GOTO 760 740 ODM=ARGS/DSQRT(1.0D0-ARGS*ARGS) DBP=DSO+JSIGN*DATAN(ODM) 760 POM=DSIN(DBP/2.0D0) DD=DSQRT(DQ*DQ+DE*DE+DAI+(EA(2)+EB(2))*(EA(2)+EB(2))*POM*POM) RETURN C END C --------------------------------------------------------------------- SUBROUTINE APPROACH(TMIN,TMAX,EA,WM1,WM2) C --------------------------------------------------------------------- C - CALCULATION OF SOME QUANTITIES CHARACTERIZING THE MEAN ORBIT OF C THE METEOR SHOWER C C INPUT: C "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE C POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 1950.0. C ORBITAL ELEMENTS OF THE ORBIT OF THE PARENT BODY 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 YEAR, MONTH, DAY - TIME OF PERIHELION PASSAGE OF THE PARENT BODY C (CONVERTED TO: EA(6) - THE TIME IN JULIAN CENTURIES FROM C JD = 2433282.423) C C OUTPUT: C SOME QUANTITIES DESCRIBING THEORETICAL RADIANTS AT BOTH ARCS (PRE- C PERIHELION AND POST-PERIHELION) OF THE ORBIT OF THE PARENT BODY; C THE OUTPUT QUANTITIES AT THE POST-PERIHELION ARC ARE: C WM1(1) - MOMENT OF MINIMUM DISTANCE BETWEEN THE EARTH AND THE C ORBIT OF THE PARENT BODY [JULIAN CENTURIES FROM 1950.0] C WM1(2) - MINIMUM DISTANCE BETWEEN POST-PERIHELION ARC OF THE C PARENT BODY ORBIT AND EARTH'S ORBIT [AU] C WM1(3) - DISTANCE OF THE PARENT BODY FROM THE EARTH AT THE MOMENT C OF MAXIMUM [AU] C WM1(4) - TIME INTERVAL BETWEEN THE PASSAGES OF THE PARENT BODY C AND EARTH ACROSS THE NEAREST POINTS OF POST-PERIHELION C ARC OF THE PARENT BODY ORBIT AND EARTH'S ORBIT [DAYS] C WM1(5) - TRUE ANOMALY OF THE PARENT BODY IN THE MOMENT OF ITS C NEAREST POST-PERIHELION APPROACH TO THE EARTH'S ORBIT C [DEGREES] C WM1(6) - TRUE ANOMALY OF THE EARTH IN THE MOMENT OF ITS NEAREST C APPROACH TO THE POST-PERIHELION ARC OF THE PARENT BODY C ORBIT [DEGREES] C THE OUTPUT QUANTITIES AT THE PRE-PERIHELION ARC ARE ANALOGOUSLY C STORED IN FIELD "WM2(I)", WHERE "I = 1, 2, 3,..., 6". C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I,J) DIMENSION EA(6),WM1(6),WM2(6),WEL(7) Cx ,EZ(6) C PI=4.0D0*DATAN(1.0D0) PI180=PI/0.18D3 GAUSS=.17202098950D-1 EPS=2.3439291D1*PI180 CEPS=DCOS(EPS) SEPS=DSIN(EPS) NTARG=3 NCENT=11 C TMN1=(TMIN+TMAX)/2.0D0 TMN2=TMN1 DT0=(TMAX-TMIN)/3.6525D2 C CHO=DCOS(PI180*EA(4)) SHO=DSIN(PI180*EA(4)) CSK=DCOS(PI180*EA(5)) SSK=DSIN(PI180*EA(5)) PRM=EA(1)*(1.0D0+EA(2)) IF(EA(2).GT.0.99999999D0) GOTO 30 A=EA(1)/(1.0D0-EA(2)) PP05=PI*A*DSQRT(A)/GAUSS PP=2.0D0*PP05 C C === 1-ST ARC === C 30 VKD=0.0D0 VKU=0.179D3 DVK=1.0D0 TA=TMIN TB=TMAX DT=DT0 C 50 WM1(2)=1.111D3 C > DO 200 VK=VKD,VKU,DVK IVKU=(VKU+1.0D-11-VKD)/DVK DO 200 IVK=0,IVKU VK=VKD+IVK*DVK WEN=1.0D0+EA(2)*DCOS(PI180*VK) RK=1.0D6 IF(DABS(WEN).LT.1.0D-12) GOTO 100 RK=PRM/WEN 100 CVS=DCOS(PI180*(VK+EA(3))) SVS=DSIN(PI180*(VK+EA(3))) XK=RK*(CVS*CHO-SVS*CSK*SHO) YK=RK*(CVS*SHO+SVS*CSK*CHO) ZK=RK*SVS*SSK C C > DO 140 ETT=TA,TB,DT IETTU=(TB+1.0D-11-TA)/DT DO 140 IETT=0,IETTU ETT=TA+IETT*DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) DR=DSQRT((XK-XZ)*(XK-XZ)+(YK-YZ)*(YK-YZ)+(ZK-ZZ)*(ZK-ZZ)) IF(DR.GT.1.0D2) GOTO 240 IF(DR.GE.WM1(2)) GOTO 140 WM1(2)=DR WM1(5)=VK Cx WM1(6)=WEL(7) C- WM1(6)=EZ(6) TMN1=ETT 140 CONTINUE 200 CONTINUE C BOTH LABELS 200 AND 240 ARE NECESSARY! 240 CONTINUE IF(DVK.GT.0.5D0) GOTO 245 IF(DABS(WM1(2)-RMOLD).LT.1.0D-5) GOTO 260 245 VKD=WM1(5)-1.0D0*DVK VKU=WM1(5)+1.0D0*DVK TA=TMN1-1.0D0*DT TB=TMN1+1.0D0*DT DVK=0.1D-1*DVK DT=0.1D-1*DT RMOLD=WM1(2) GOTO 50 C 260 CONTINUE C IF(WM1(2).GT.1.0D2) GOTO 530 Cx DTM=3.6525D4*(TMN1-EA(6)) C? DTM=TMN1-EA(6) C? READ(*,*) JXXX C? EQX=EA(1) C? EEX=EA(2) C? CALL KEPLER(EQX,EEX,DTM,QWV) C? WENK=1.0D0+EA(2)*DCOS(QWV) C? RK=PRM/WENK C? CVS=DCOS(QWV+PI180*EA(3)) C? SVS=DSIN(QWV+PI180*EA(3)) C? XK=RK*(CVS*CHO-SVS*CSK*SHO) C? YK=RK*(CVS*SHO+SVS*CSK*CHO) C? ZK=RK*SVS*SSK C? CALL PLEPH(TMN1,NTARG,NCENT,RRD) C CALL ECLXYZ(TMN1,XZ,YZ,ZZ,WEL) C? XZ=RRD(1) C? YZ=RRD(2)*CEPS+RRD(3)*SEPS C? ZZ=-RRD(2)*SEPS+RRD(3)*CEPS C VXZ=RRD(4) C VYZ=RRD(5)*CEPS+RRD(6)*SEPS C VZZ=-RRD(5)*SEPS+RRD(6)*CEPS C CALL ELEMS(ETT,XZ,YZ,VZ,VXZ,VYZ,VZZ,EZ) C? DR=DSQRT((XK-XZ)*(XK-XZ)+(YK-YZ)*(YK-YZ)+(ZK-ZZ)*(ZK-ZZ)) C? WM1(3)=DR C C??? WM1(4)=DTM C? WM1(4)=TMN1 C? QWV=WM1(5)*PI180 C? CALL TIMEANOM(EQX,EEX,QWV,DTT) C? WM1(4)=WM1(4)-DTT C? IF(EA(2).GT.0.99999999D0) GOTO 270 C? 263 CONTINUE C? IF(WM1(4).LE.PP05) GOTO 265 C? WM1(4)=WM1(4)-PP C? GOTO 263 C? 265 CONTINUE C? IF(WM1(4).GE.-PP05) GOTO 270 C? WM1(4)=WM1(4)+PP C? GOTO 265 C 270 WM1(1)=TMN1 Ct WRITE(*,*) 'TMN1, EA(6) = ',TMN1,EA(6) C C === 2-ND ARC === C 530 VKD=0.359D3 VKU=0.181D3 DVK=-1.0D0 TA=TMIN TB=TMAX DT=DT0 C 550 WM2(2)=1.111D3 C > DO 700 VK=VKD,VKU,DVK IVKU=(VKU+1.0D-11-VKD)/DVK DO 700 IVK=0,IVKU VK=VKD+IVK*DVK WEN=1.0D0+EA(2)*DCOS(PI180*VK) RK=1.0D6 IF(DABS(WEN).LT.1.0D-12) GOTO 600 RK=PRM/WEN 600 CVS=DCOS(PI180*(VK+EA(3))) SVS=DSIN(PI180*(VK+EA(3))) XK=RK*(CVS*CHO-SVS*CSK*SHO) YK=RK*(CVS*SHO+SVS*CSK*CHO) ZK=RK*SVS*SSK C C > DO 640 ETT=TA,TB,DT IETTU=(TB+1.0D-11-TA)/DT DO 640 IETT=0,IETTU ETT=TA+IETT*DT c CALL PLEPH(ETT,NTARG,NCENT,RRD) Ct WRITE(*,*) 'APPROACH-3: Von z PLEPH.' T20=(ETT-2.451545D6)/3.6525D4 CALL EARTHrv(T20,XZ,YZ,ZZ,VXZ,VYZ,VZZ,WEL) DR=DSQRT((XK-XZ)*(XK-XZ)+(YK-YZ)*(YK-YZ)+(ZK-ZZ)*(ZK-ZZ)) IF(DR.GT.1.0D2) GOTO 740 IF(DR.GE.WM2(2)) GOTO 640 WM2(2)=DR WM2(5)=VK Cx WM2(6)=WEL(7) C- WM2(6)=EZ(6) TMN2=ETT 640 CONTINUE 700 CONTINUE C BOTH LABELS 700 AND 740 ARE NECESSARY! 740 CONTINUE IF(DVK.GT.0.5D0) GOTO 745 IF(DABS(WM2(2)-RMOLD).LT.1.0D-5) GOTO 760 745 VKD=WM2(5)-1.0D0*DVK VKU=WM2(5)+1.0D0*DVK TA=TMN2-1.0D0*DT TB=TMN2+1.0D0*DT DVK=0.1D-1*DVK DT=0.1D-1*DT RMOLD=WM2(2) GOTO 550 C 760 CONTINUE C IF(WM2(2).GT.1.0D2) GOTO 800 Cx DTM=3.6525D4*(TMN2-EA(6)) C? DTM=TMN2-EA(6) C? EQX=EA(1) C? EEX=EA(2) C? CALL KEPLER(EQX,EEX,DTM,QWV) C? WENK=1.0D0+EA(2)*DCOS(QWV) C? RK=PRM/WENK C? CVS=DCOS(QWV+PI180*EA(3)) C? SVS=DSIN(QWV+PI180*EA(3)) C? XK=RK*(CVS*CHO-SVS*CSK*SHO) C? YK=RK*(CVS*SHO+SVS*CSK*CHO) C? ZK=RK*SVS*SSK C? CALL PLEPH(TMN2,NTARG,NCENT,RRD) C CALL ECLXYZ(TMN2,XZ,YZ,ZZ,WEL) C? XZ=RRD(1) C? YZ=RRD(2)*CEPS+RRD(3)*SEPS C? ZZ=-RRD(2)*SEPS+RRD(3)*CEPS C VXZ=RRD(4) C VYZ=RRD(5)*CEPS+RRD(6)*SEPS C VZZ=-RRD(5)*SEPS+RRD(6)*CEPS C CALL ELEMS(ETT,XZ,YZ,VZ,VXZ,VYZ,VZZ,EZ) C? DR=DSQRT((XK-XZ)*(XK-XZ)+(YK-YZ)*(YK-YZ)+(ZK-ZZ)*(ZK-ZZ)) C? WM2(3)=DR C C??? WM2(4)=DTM C? WM2(4)=TMN2 C? QWV=WM2(5)*PI180 C? CALL TIMEANOM(EQX,EEX,QWV,DTT) C? WM2(4)=WM2(4)-DTT C? IF(EA(2).GT.0.99999999D0) GOTO 770 C? 763 CONTINUE C? IF(WM2(4).LE.PP05) GOTO 765 C? WM2(4)=WM2(4)-PP C? GOTO 763 C? 765 CONTINUE C? IF(WM2(4).GE.-PP05) GOTO 770 C? WM2(4)=WM2(4)+PP C? GOTO 765 C 770 WM2(1)=TMN2 800 CONTINUE RETURN C END C -------------------------------------------------------------------- SUBROUTINE JULIAN(IROK,IMES,DEN,T20,ETT) C -------------------------------------------------------------------- C - CALCULATION OF JULIAN DATE FROM CIVIL DATE C C INPUT: C IROK, IMES, DEN - YEAR, MONTH, AND DAY OF THE DATE C C OUTPUT: C T20 - GIVEN TIME ACCOUNTED FROM THE BEGINNING OF EPOCH 1900.0 C (JD = 2415020.0) AND EXPRESSED IN JULIAN CENTURIES C ETT - JDT OF THE GIVEN DATE C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8 (A-H,K-Z) IMPLICIT INTEGER (I,J) DIMENSION JMESI(13) DATA JMESI/0,31,28,31,30,31,30,31,31,30,31,30,31/ C IF(IROK.GT.1582) GOTO 101 IF(IROK.EQ.1582) GOTO 102 IF(IROK.GT.0.AND.IROK.LT.1582) GOTO 103 IF(IROK.GE.-4713.AND.IROK.LT.0) GOTO 104 C ETT=(IROK+4713)*365 ETT=ETT-INT((-IROK-4712)/4) 100 T0JD=-0.5D0 GOTO 210 101 ETT=(IROK-1583)*365 ETT=ETT+INT((IROK-1581)/4) ETT=ETT-INT((IROK-1501)/100)+INT((IROK-1201)/400) T0JD=2.2992385D6 GOTO 210 102 ETT=0.0D0 IF(IMES.EQ.10.AND.DEN.GT.1.49999999D1) ETT=-1.0D1 IF(IMES.GT.10) ETT=-1.0D1 T0JD=2.2988835D6 GOTO 210 103 ETT=(IROK-1)*365 ETT=ETT+INT((IROK-1)/4) T0JD=1.7214235D6 GOTO 210 104 ETT=(IROK+4713)*365 ETT=ETT+INT((IROK+4715)/4) T0JD=-0.5D0 210 YY=IROK/1.00D2-INT(IROK/100) ZZ=IROK/4.00D2-INT(IROK/400) IF(IROK.GT.1582.AND.YY.EQ.0.0D0.AND.ZZ.GT.0.0D0) GOTO 220 YY=IROK/4.0D0-INT(IROK/4) IF(YY.EQ.0.0D0) JMESI(3)=29 IF(IROK.EQ.-1) JMESI(3)=29 220 CONTINUE DO 225 JJM=1,IMES-1 ETT=ETT+JMESI(JJM+1) 225 CONTINUE ETT=ETT+DEN-1.0D0+T0JD JMESI(3)=28 C T20=(ETT-2.415020D6)/3.6525D4 CC ETT=(ETT-2.433282423D6)/3.6525D4 RETURN C END C --------------------------------------------------------------------- SUBROUTINE INVJUL(TJD,IRK,IMS,DNI) C --------------------------------------------------------------------- C - CALCULATION OF CIVIL DATE FROM JULIAN DATE (INVERSION OPERATION C TO THAT EXECUTED BY SUBROUTINE "JULIAN") C C INPUT: C TJD - TIME FROM THE BEGINNING OF EPOCH 2000.0 (JD = 2433282.423) C IN JULIAN CENTURIES (IT CAN BE CALCULATED FROM JULIAN DATE C "JD" IN DAYS AS: TJD=(JD-2433282.423)/36525) C C OUTPUT: C IRK, IMS, DNI - YEAR, MONTH, AND DAY OF THE CORRESPONDING CIVIL C DATE C C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C IMPLICIT REAL*8 (D,T) IMPLICIT INTEGER (I) DIMENSION IMESI(13) DATA IMESI/0,31,28,31,30,31,30,31,31,30,31,30,31/ C CC DNI=TJD*3.6525D4+2.433282423D6+1.0D0 DNI=TJD+1.0D0 IF(DNI.LT.0.5D0) GOTO 16 IF(DNI.GE.0.5D0.AND.DNI.LT.1.7214245D6) GOTO 15 IF(DNI.GE.1.7214245D6.AND.DNI.LT.2.2991615D6) GOTO 14 IF(DNI.GE.2.2991615D6.AND.DNI.LT.2.2992395D6) GOTO 13 IF(DNI.GE.2.2992395D6.AND.DNI.LT.2.4515455D6) GOTO 12 C C FROM 2000-JAN-1.0 11 CONTINUE DNI=DNI-2.4515445D6 IRK=INT(DNI/3.65242198D2) DN0=365*IRK+INT((IRK+3)/4.0D0)-INT((IRK-1)/1.0D2)+INT((IRK-1)/4.0D *2) DNI=DNI-DN0 IF(DNI.LT.1.0D0) GOTO 30 IRK=IRK+2000 GOTO 20 30 IRK=IRK-1 DNI=DNI+3.65D2 DH4=IRK/4.0D0-INT(IRK/4.0D0) DH100=IRK/1.0D2-INT(IRK/1.0D2) IF(DABS(DH4).LT.1.0D-6.AND.DH100.GT.0.0D0) DNI=DNI+1.0D0 DH400=IRK/4.0D2-INT(IRK/4.0D2) IF(DABS(DH400).LT.1.0D-6) DNI=DNI+1.0D0 IRK=IRK+2000 GOTO 20 C C FROM 1583-JAN-1.0 TO 2000-JAN-1.0 12 CONTINUE DNI=DNI-2.2992385D6 IRK=INT(DNI/3.65242198D2) DN0=365*IRK+INT((IRK+2)/4.0D0)-INT((IRK+82)/1.0D2)+INT((IRK+382)/4 *.0D2) DNI=DNI-DN0 IF(DNI.LT.1.0D0) GOTO 35 IRK=IRK+1583 GOTO 20 35 IRK=IRK-1 DH4=(IRK+3)/4.0D0-INT((IRK+3)/4.0D0) DH100=(IRK+83)/1.0D2-INT((IRK+83)/1.0D2) DNI=DNI+3.65D2 IF(DABS(DH4).LT.1.0D-6.AND.DH100.GT.0.0D0) DNI=DNI+1.0D0 DH400=(IRK+383)/4.0D2-INT((IRK+383)/4.0D2) IF(DABS(DH400).LT.1.0D-6) DNI=DNI+1.0D0 IRK=IRK+1583 GOTO 20 C C FROM 1582-OCT-15.0 TO 1583-JAN-1.0 13 CONTINUE IRK=1582 DNI=DNI-2.2991605D6+1.4D1 IF(DNI.GE.3.2D1) GOTO 40 IMS=10 GOTO 26 40 CONTINUE IF(DNI.GE.6.2D1) GOTO 45 DNI=DNI-3.1D1 IMS=11 GOTO 26 45 DNI=DNI-6.1D1 IMS=12 GOTO 26 C C FROM 1-JAN-1.0 TO 1582-OCT-15.0 14 CONTINUE DNI=DNI-1.7214235D6 IRK=INT(DNI/3.6525D2)+1 DN0=3.65D2*(IRK-1)+INT((IRK-1)/4.0D0) DNI=DNI-DN0 IF(DNI.GE.1.0D0) GOTO 20 IRK=IRK-1 DH4=IRK/4.0D0-INT(IRK/4.0D0) DNI=DNI+3.65D2 IF(DABS(DH4).LT.1.0D-6) DNI=DNI+1.0D0 GOTO 20 C C FROM -4713-JAN-1.0 TO 1-JAN-1.0 15 CONTINUE DNI=DNI+0.5D0 IRK=INT(DNI/3.6525D2) DN0=365*IRK+INT((IRK+2)/4.0D0) IF(DNI.GT.1.721118D6) DNI=DNI-1.0D0 DNI=DNI-DN0 IF(DNI.LT.1.0D0) GOTO 50 IRK=IRK-4713 GOTO 20 50 IRK=IRK-1 DH4=(IRK+3)/4.0D0-INT((IRK+3)/4.0D0) DNI=DNI+3.65D2 IF(DABS(DH4).LT.1.0D-6) DNI=DNI+1.0D0 IRK=IRK-4713 GOTO 20 C C TO -4713-JAN-1.0 16 CONTINUE DNI=DNI+1.9307115D6 IRK=INT(DNI/3.6525D2) DN0=365*IRK+INT(IRK/4.0D0) DNI=DNI-DN0 IF(DNI.LT.1.0D0) GOTO 55 IRK=IRK-9999 GOTO 20 55 IRK=IRK-1 DH4=(IRK+1)/4.0D0-INT((IRK+1)/4.0D0) DNI=DNI+3.65D2 IF(DABS(DH4).LT.1.0D-6) DNI=DNI+1.0D0 IRK=IRK-9999 C 20 CONTINUE DH4=IRK/4.0D0-INT(IRK/4.0D0) DH100=IRK/1.0D2-INT(IRK/1.0D2) DH400=IRK/4.0D2-INT(IRK/4.0D2) IF(DABS(DH4).LT.1.0D-6) IMESI(3)=29 IF(IRK.GT.1582.AND.DABS(DH100).LT.1.0D-6) IMESI(3)=28 IF(IRK.GT.1582.AND.DABS(DH400).LT.1.0D-6) IMESI(3)=29 IF(IRK.EQ.0) IRK=-1 IMS=1 23 CONTINUE IF(IMS.NE.13) GOTO 60 IMS=12 DNI=DNI-1.0D0 60 CONTINUE IF(IMESI(IMS+1)+1.GT.DNI) GOTO 26 DNI=DNI-IMESI(IMS+1) IMS=IMS+1 GOTO 23 26 CONTINUE IMESI(3)=28 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