C23456789012345678901234567890123456789012345678901234567890123456789012
C LAST MODIFIED SEPTEMBER 16, 1988
C
C *** BINGHAM ***
C
C SUBROUTINE TO GENERATE POINTS ON THE CONFIDENCE REGION BOUNDARY, FOR
C USE BY SURFACE2, OR OTHER CONTOURING PROGRAM.
C INPUT:
C Q - THE MATRIX Q (IN SYMMETRIC STORAGE MODE)
C LHAT - THE SMALLEST EIGENVALUE OF Q
C FK - CONSTANT USED IN DETERMINING THE CONFIDENCE REGION
C UXX - OPTIMAL AXIS OF ROTATION
C OUTPUT:
C JER - ERROR INDICATOR
C 0 - ALL IS WELL
C 1 - AN ERROR HAS OCCURRED
C OUTPUT FILES:
C 30 - BOUNDING CURVE
C 31 - UPPER SURFACE
C 32 - LOWER SURFACE
SUBROUTINE BINGHAM(Q,LHAT,FK,UXX,JER)
REAL Q(10),QF(4,4),NU(3),W(3,3),QT(10)
REAL UXX(3),LHAT,MF(3,3)
CHARACTER SW*1
DATA PI/3.1415927/
PRINT *,' '
PRINT *,'DETERMINING THE CONFIDENCE REGION IN THE FORM: MIN. AND'
PRINT *,'MAX. ANGLES OF ROTATION EXPRESSED AS FUNCTIONS OF AXIS'
PRINT *,'LONGITUDE AND AXIS LATITUDE.'
PRINT *,' '
IND=0
PRINT *,'DO YOU WANT OUTPUT FILES PRODUCED SUITABLE FOR GRAPHING'
PRINT *,'THE CONFIDENCE REGION?'
READ 999,SW
IF ((SW.EQ.'Y').OR.(SW.EQ.'y')) IND=1
JER=0
C MODIFY THE MATRIX Q BY SUBTRACTING LHAT TIMES THE IDENTITY FROM Q,
C AND DIVIDING THE RESULT BY FK.
DO 10 I=1,10
10 QT(I)=Q(I)
QT(1)=QT(1) - LHAT
QT(3)=QT(3) - LHAT
QT(6)=QT(6) - LHAT
QT(10)=QT(10) - LHAT
DO 20 I=1,10
20 QT(I)=QT(I)/FK
K=0
DO 32 I=1,4
DO 32 J=1,I
K=K+1
QF(I,J)=QT(K)
32 QF(J,I)=QF(I,J)
PRINT *,'THE MATRIX QTILDE'
DO 35 I=1,4
35 PRINT 990,(QF(I,J),J=1,4)
990 FORMAT (' ',4E18.8)
C TEST WHETHER THE IDENTITY IS IN THE CONFIDENCE REGION. IF NOT,
C CALCULATE THE MATRIX
C M = (AZERO - 1)*B - DVECT*TRANSPOSE(DVECT)
C AND FIND ITS EIGENVALUES AND EIGENVECTORS.
AZERO=QF(1,1)
IF (AZERO.LE.1.0) THEN
ICASE=7
ELSE
CALL MEIG(QF,UXX,NU,W,MF,ICASE,IER)
IF (IER.NE.0) THEN
JER=1
RETURN
END IF
END IF
C PRINT A DESCRIPTION OF THE SET A OF ADMISSIBLE AXES, DEPENDING ON
C THE VALUE OF ICASE.
IF (ICASE.EQ.1) THEN
PRINT *,' '
PRINT *,'THE SET A OF ADMISSIBLE AXES IS A CAP NOT CONTAINING'
PRINT *,'EITHER POLE. IN THE AXIS LONGITUDE-AXIS LATITUDE'
PRINT *,'PLANE THIS SET IS BOUNDED BY A CLOSED CURVE.'
ELSE IF (ICASE.EQ.2) THEN
PRINT *,' '
PRINT *,'THE SET A OF ADMISSIBLE AXES IS A CAP CONTAINING'
PRINT *,'THE NORTH POLE. IN THE AXIS LONGITUDE-AXIS LATITUDE'
PRINT *,'PLANE THIS SET IS BOUNDED BY THE LINES: AXIS LONGI-'
PRINT *,'TUDE = -180 DEGREES, AXIS LONGITUDE = 180 DEGREES,'
PRINT *,'AXIS LATITUDE = 90 DEGREES, AND A CURVE WHICH FORMS'
PRINT *,'THE SOUTHERN BORDER.'
ELSE IF (ICASE.EQ.3) THEN
PRINT *,' '
PRINT *,'THE SET A OF ADMISSIBLE AXES IS A CAP CONTAINING'
PRINT *,'THE SOUTH POLE. IN THE AXIS LONGITUDE-AXIS LATITUDE'
PRINT *,'PLANE THIS SET IS BOUNDED BY THE LINES: AXIS LONGI-'
PRINT *,'TUDE = -180 DEGREES, AXIS LONGITUDE = 180 DEGREES,'
PRINT *,'AXIS LATITUDE = -90 DEGREES, AND A CURVE WHICH FORMS'
PRINT *,'THE NORTHERN BORDER.'
ELSE IF (ICASE.EQ.4) THEN
PRINT *,' '
PRINT *,'THE SET A OF ADMISSIBLE AXES IS THE COMPLEMENT OF'
PRINT *,'TWO ANTI-PODAL CAPS WHICH CONTAIN THE POLES.'
PRINT *,'HENCE, THIS SET IS AN EQUATORIAL BELT. IN THE AXIS'
PRINT *,'LONGITUDE-AXIS LATITUDE PLANE THIS SET IS BOUNDED BY'
PRINT *,'THE LINES: AXIS LONGITUDE = -180 DEGREES, AXIS LONGI-'
PRINT *,'TUDE = 180 DEGREES, AND TWO CURVES WHICH FORM THE'
PRINT *,'NORTHERN AND SOUTHERN BORDERS OF THE BELT.'
ELSE IF (ICASE.EQ.5) THEN
PRINT *,' '
PRINT *,'THE SET A OF ADMISSIBLE AXES IS THE COMPLEMENT OF'
PRINT *,'TWO ANTI-PODAL CAPS WHICH DO NOT CONTAIN THE'
PRINT *,'POLES. IN THE AXIS LONGITUDE-AXIS LATITUDE PLANE'
PRINT *,'THIS SET IS BOUNDED BY THE FOUR LINES: AXIS LONGI-'
PRINT *,'TUDE = -180 DEGREES, AXIS LONGITUDE = 180 DEGREES,'
PRINT *,'AXIS LATITUDE = -90 DEGREES, AND AXIS LATITUDE ='
PRINT *,'90 DEGREES; HOWEVER, THERE ARE TWO HOLES IN THE SET.'
ELSE IF (ICASE.EQ.6) THEN
PRINT *,' '
PRINT *,'ANY AXIS IS ADMISSIBLE; HOWEVER, THE IDENTITY IS NOT'
PRINT *,'IN THE CONFIDENCE REGION. THE SET A OF ADMISSIBLE'
PRINT *,'AXES IS THE ENTIRE AXIS LONGITUDE-AXIS LATITUDE PLANE.'
ELSE IF (ICASE.EQ.7) THEN
PRINT *,' '
PRINT *,'THE IDENTITY IS IN THE CONFIDENCE REGION. HENCE, EACH'
PRINT *,'AXIS IS ADMISSIBLE. THAT IS, THE SET A OF ADMISSIBLE'
PRINT *,'AXES IS THE ENTIRE AXIS LONGITUDE-AXIS LATITUDE PLANE.'
END IF
CALL BOUNDC(AZERO,NU,W,ICASE,IND,PMIND,PMAXD,TMIND,TMAXD,IER)
IF (IER.NE.0) THEN
JER=1
RETURN
END IF
PRINT *,' '
PRINT *,'TYPE C TO CONTINUE.'
READ 999,SW
999 FORMAT (A)
ALAT1=FLOAT(LARGIN(TMIND))
ALAT2=FLOAT(-LARGIN(-TMAXD))
ALONG1=FLOAT(LARGIN(PMIND))
ALONG2=FLOAT(-LARGIN(-PMAXD))
CALL GRID(ALAT1,ALAT2,ALONG1,ALONG2,QF,MF,W,ICASE,IND,IER)
IF (IER.NE.0) THEN
JER=1
RETURN
END IF
RETURN
END
C
C *** BOUNDC ***
C
C SUBROUTINE TO CALCULATE POINTS ON THE BOUNDING CURVE.
C INPUTS:
C AZERO - QF(1,1)
C NU - ARRAY CONTAINING THE EIGENVALUES OF M
C W - ARRAY CONTAINING THE EIGENVECTORS OF M
C ICASE - INDICATOR OF THE TYPE OF SET A
C IND - INDICATOR
C IND=0: BOX OUT THE CONFIDENCE REGION
C IND=1: WRITE THE BOUNDARY CONTOUR TO FILE 30
C OUTPUTS:
C PMIND, PMAXD - MIN. AND MAX. LONGITUDES OVER THE CONFIDENCE REGION
C TMIND, TMAXD - MIN. AND MAX. LATITUDES OVER THE CONFIDENCE REGION
C JER - ERROR INDICATOR
C 0 - ALL IS WELL
C 1 - AN ERROR HAS OCCURRED
C DEFINED BY PARAMETER STATEMENTS:
C NPART - NO. OF PARTS INTO WHICH TO DIVIDE THE BOUNDING CURVE,
C TO GET THE STEP-SIZE
C KMAX - MAX. NO. OF POINTS (MUST BE AT LEAST 13)
C THE SUBROUTINE BOUNDC WRITES THE LONGITUDE AND LATITUDE VALUES FOR THE
C POINTS ON THE BOUNDING CURVE TO FILE 30, PROVIDED IND=1 AND
C ICASE.LE.5.
SUBROUTINE BOUNDC(AZERO,NU,W,ICASE,IND,PMIND,PMAXD,TMIND,TMAXD,
& JER)
PARAMETER (NPART=300)
PARAMETER (KMAX=400)
REAL NU(3),W(3,3)
REAL PHI(KMAX),THETA(KMAX),PHIT(KMAX)
REAL PHIM(KMAX),THETM(KMAX)
REAL U(3),OLDU(3)
CHARACTER BFILE*10
DATA PI/3.1415927/
JER=0
IF (ICASE.GE.6) THEN
PMIND= -180.0
PMAXD= 180.0
TMIND= -90.0
TMAXD= 90.0
PRINT 811,PMIND,PMAXD,TMIND,TMAXD
RETURN
END IF
IF (IND.EQ.1) THEN
PRINT *,' TYPE NAME OF BOUNDARY FILE, I.E., FILE TO RECEIVE'
PRINT *,' POINTS ON THE BOUNDARY OF THE SET OF ADMISSIBLE'
PRINT *,' AXES. (MAX.=10 CHARACTERS.)'
READ 999,BFILE
999 FORMAT (A)
OPEN (30,FILE=BFILE,STATUS='UNKNOWN')
REWIND 30
END IF
C CALCULATION OF ESTIMATED LENGTH OF THE BOUNDING CURVE. THIS IS
C DONE BY CALCULATING THIRTEEN POINTS AROUND THE CURVE, AT EQUAL
C INCREMENTS IN PHIT.
DELPT=PI/6.0
DO 105 K=1,13
PHIT(K)=(K-1)*DELPT
CALL EVALF(PHIT(K),NU,W,AZERO,FVAL,PHI(K),THETA(K),U,K-1,
&OLDU,OPHI,IER)
IF (IER.NE.0) THEN
JER=1
RETURN
END IF
DO 104 I=1,3
104 OLDU(I)=U(I)
105 OPHI=PHI(K)
ELENG=0.0
DO 107 K=1,12
DIST=SQRT((PHI(K+1) - PHI(K))*(PHI(K+1) - PHI(K))
& + (THETA(K+1) - THETA(K))*(THETA(K+1) - THETA(K)))
107 ELENG=ELENG + DIST
ELENGD=ELENG*180./PI
DELS=ELENG/FLOAT(NPART)
PHIT(1)=0.0
IF (ICASE.EQ.3) DELS= -DELS
C MAIN LOOP. CALCULATE THE DESIRED PHIT VALUES BY APPLYING EULER'S
C METHOD TO THE DIFFERENTIAL EQUATION:
C DPHIT/DS = F(PHIT) ,
C WHERE F(PHIT) = 1./SQRT((DPHI/DPHIT)**2 + (DTHETA/DPHIT)**2).
C HERE S IS THE ARC LENGTH ALONG THE BOUNDING CURVE, MEASURED IN
C RADIANS.
DO 120 K=2,KMAX
CALL EVALF(PHIT(K-1),NU,W,AZERO,FVAL,PHI(K-1),THETA(K-1),U,
&K-2,OLDU,OPHI,IER)
IF (IER.NE.0) THEN
JER=1
RETURN
END IF
DO 118 I=1,3
118 OLDU(I)=U(I)
OPHI=PHI(K-1)
PHIT(K)=PHIT(K-1) + FVAL*DELS
IF (ABS(PHIT(K)).GE.(2.*PI)) GO TO 122
120 CONTINUE
C END OF MAIN LOOP.
JER=1
PRINT *,'ERROR IN BOUNDC:'
PRINT *,'UNABLE TO CLOSE THE BOUNDING CURVE WITH THE PRESCRIBED'
PRINT *,'NPART AND KMAX. EITHER NPART NEEDS TO BE DECREASED OR'
PRINT *,'KMAX NEEDS TO BE INCREASED.'
RETURN
C THE BOUNDING CURVE IS CLOSED. BRANCH TO DIFFERENT PROCEDURES
C DEPENDING ON THE VALUE OF ICASE.
122 CONTINUE
IF ((ICASE.GE.2).AND.(ICASE.LE.4)) GO TO 200
IF (ICASE.EQ.5) GO TO 500
C ICASE = 1. THE SET A IS A CAP NOT CONTAINING EITHER POLE.
C CENTER THE PHI VALUES, IF NECESSARY.
NPT=K
PHIT(NPT)=2.*PI
PHI(NPT)=PHI(1)
THETA(NPT)=THETA(1)
PMIN=PHI(1)
PMAX=PHI(1)
DO 125 K=1,NPT
IF (PHI(K).LT.PMIN) PMIN=PHI(K)
IF (PHI(K).GT.PMAX) PMAX=PHI(K)
125 CONTINUE
PMEAN=(PMIN + PMAX)/2.
IF (PMEAN.GT.PI) THEN
DO 128 K=1,NPT
128 PHI(K)=PHI(K) - 2.*PI
PMIN=PMIN - 2.*PI
PMAX=PMAX - 2.*PI
ELSE IF (PMEAN.LE.(-PI)) THEN
DO 133 K=1,NPT
133 PHI(K)=PHI(K) + 2.*PI
PMIN=PMIN + 2.*PI
PMAX=PMAX + 2.*PI
END IF
TMIN=THETA(1)
TMAX=THETA(1)
DO 172 K=1,NPT
IF (THETA(K).LT.TMIN) TMIN=THETA(K)
IF (THETA(K).GT.TMAX) TMAX=THETA(K)
172 CONTINUE
PMIND=PMIN*180./PI
PMAXD=PMAX*180./PI
TMIND=TMIN*180./PI
TMAXD=TMAX*180./PI
PRINT 811,PMIND,PMAXD,TMIND,TMAXD
811 FORMAT (' ','MIN., MAX. AXIS LONGITUDE OVER THE SET A:'/' ',5X,
&F10.5,',',F10.5,' DEGREES'/' ','MIN., MAX. AXIS LATITUDE OVER',
&' THE SET A:'/' ',5X,F10.5,',',F10.5,' DEGREES')
IF (IND.EQ.1) THEN
WRITE (30,801) NPT
801 FORMAT (I4,' 1')
DO 180 K=1,NPT
BLONG=PHI(K)*180./PI
BLAT=THETA(K)*180./PI
180 WRITE (30,802) BLONG,BLAT
802 FORMAT (1X,2E16.8)
CLOSE (30)
PRINT 703,NPT,BFILE
703 FORMAT (' ',' A SEQUENCE OF ',I3,' POINTS AROUND THE CLOSED',
&' CURVE BOUNDING'/' ',' THE SET OF ADMISSIBLE AXES HAS',
&' BEEN WRITTEN TO'/' ',' THE FILE ',A10,'.')
END IF
RETURN
C ICASE = 2, 3, OR 4.
C ICASE = 2: THE SET A IS A CAP CONTAINING THE NORTH POLE.
C ICASE = 3: THE SET A IS A CAP CONTAINING THE SOUTH POLE.
C ICASE = 4: THE SET A IS THE COMPLEMENT OF TWO ANTI-PODAL CAPS
C WHICH CONTAIN THE POLES.
C CALCULATE A SEQUENCE OF POINTS (PHIM(K),THETM(K),K=1,2,...,NPT,
C WHICH ARE ARRANGED IN ORDER ON PHI, WITH PHI GOING FROM -180 TO
C 180.
200 CONTINUE
NPT=K - 1
DO 205 K=2,NPT
IF (PHI(K).GT.PI) GO TO 215
205 CONTINUE
DO 210 K=1,NPT
PHIM(K)=PHI(K)
210 THETM(K)=THETA(K)
GO TO 225
215 CONTINUE
KZ=K-1
DO 220 L=1,NPT-KZ
PHIM(L)=PHI(KZ + L) - 2.*PI
220 THETM(L)=THETA(KZ + L)
DO 222 L=1,KZ
K=L + NPT - KZ
PHIM(K)=PHI(L)
222 THETM(K)=THETA(L)
C FROM NOW ON WE TREAT THE CASES ICASE=2,3 AND ICASE=4 DIFFERENTLY.
225 CONTINUE
IF (ICASE.EQ.4) GO TO 400
TMIN=THETM(1)
TMAX=THETM(1)
DO 245 K=1,NPT
IF (THETM(K).LT.TMIN) TMIN=THETM(K)
IF (THETM(K).GT.TMAX) TMAX=THETM(K)
245 CONTINUE
PMIND= -180.
PMAXD= 180.
IF (ICASE.EQ.2) THEN
TMIND=TMIN*180./PI
TMAXD=90.0
ELSE IF (ICASE.EQ.3) THEN
TMIND= -90.0
TMAXD=TMAX*180./PI
END IF
PRINT 811,PMIND,PMAXD,TMIND,TMAXD
IF (IND.EQ.1) THEN
WRITE (30,803) NPT
803 FORMAT (I4,' 0')
DO 250 K=1,NPT
BLONG=PHIM(K)*180./PI
BLAT=THETM(K)*180./PI
250 WRITE (30,802) BLONG,BLAT
CLOSE (30)
IF (ICASE.EQ.2) THEN
PRINT 712,NPT,BFILE
ELSE IF (ICASE.EQ.3) THEN
PRINT 713,NPT,BFILE
END IF
END IF
712 FORMAT (' ',' A SEQUENCE OF ',I3,' POINTS ALONG THE SOUTHERN',
&' BORDER'/' ',' OF THE SET OF ADMISSIBLE AXES HAS BEEN',
&' WRITTEN TO'/' ',' THE FILE ',A10,'.')
713 FORMAT (' ',' A SEQUENCE OF ',I3,' POINTS ALONG THE NORTHERN',
&' BORDER'/' ',' OF THE SET OF ADMISSIBLE AXES HAS BEEN',
&' WRITTEN TO'/' ',' THE FILE ',A10,'.')
RETURN
C ICASE = 4. THE SET A IS THE COMPLEMENT OF TWO ANTI-PODAL CAPS
C WHICH CONTAIN THE POLES. KZZ IS THE NUMBER OF POINTS (PHIM(K),
C THETM(K)) WITH PHI NEGATIVE.
400 CONTINUE
DO 405 K=1,NPT
IF (PHIM(K).GT.0.) GO TO 410
405 CONTINUE
JER=1
PRINT *,'ERROR IN BOUNDC (ICASE=4):'
PRINT *,'THERE ARE NO POINTS ON THE BOUNDING CURVE WITH POSITVE'
PRINT *,'LONGITUDE. THE POINTS ON THE BOUNDING CURVE ARE TOO'
PRINT *,'SPARSE.'
RETURN
410 CONTINUE
KZZ=K-1
IF (KZZ.LE.0) THEN
JER=1
PRINT *,'ERROR IN BOUNDC (ICASE=4):'
PRINT *,'THERE ARE NO POINTS ON THE BOUNDING CURVE WITH'
PRINT *,'NEGATIVE LONGITUDE. THE POINTS ON THE BOUNDING'
PRINT *,'CURVE ARE TOO SPARSE.'
RETURN
END IF
TMAX=THETM(1)
DO 450 K=1,NPT
IF (THETM(K).GT.TMAX) TMAX=THETM(K)
450 CONTINUE
PMIND= -180.
PMAXD= 180.
TMAXD=TMAX*180./PI
TMIND= -TMAXD
PRINT 811,PMIND,PMAXD,TMIND,TMAXD
IF (IND.EQ.1) THEN
WRITE (30,803) NPT
DO 460 K=1,NPT
BLONG=PHIM(K)*180./PI
BLAT=THETM(K)*180./PI
460 WRITE (30,802) BLONG,BLAT
WRITE (30,803) NPT
DO 465 L=KZZ+1,NPT
BLONG=PHIM(L)*180./PI - 180.
BLAT= -THETM(L)*180./PI
465 WRITE (30,802) BLONG,BLAT
DO 468 L=1,KZZ
BLONG=PHIM(L)*180./PI + 180.
BLAT= -THETM(L)*180./PI
468 WRITE (30,802) BLONG,BLAT
CLOSE (30)
PRINT 714,NPT,BFILE
END IF
714 FORMAT (' ',' A SEQUENCE OF ',I3,' POINTS ALONG THE CURVE',
&' WHICH FORMS THE'/' ',' NORTHERN BORDER OF THE BELT HAVE BEEN',
&' WRITTEN TO THE FILE ',A10,'.'/' ',' THESE ARE FOLLOWED BY THE',
&' SAME NUMBER OF POINTS ALONG THE'/' ',' CURVE WHICH FORMS THE',
&' SOUTHERN BORDER.')
RETURN
C ICASE = 5. THE SET A IS THE COMPLEMENT OF TWO ANTI-PODAL CAPS
C WHICH DO NOT CONTAIN THE POLES. IF THE CAP WHICH WE HAVE CROSSES
C EITHER OF THE LINES PHI=PI OR PHI= -PI, WE REPLACE IT BY THE
C ANTI-PODAL CAP. WE THEN HAVE SUBCASE A OR SUBCASE B, DEPENDING
C ON WHETHER THE CAP WHICH WE HAVE OBTAINED DOES NOT CROSS THE
C LINE PHI=0, OR DOES CROSS IT.
500 CONTINUE
NPT=K
PHI(NPT)=PHI(1)
THETA(NPT)=THETA(1)
PMIND= -180.
PMAXD= 180.
TMIND= -90.
TMAXD= 90.
PRINT 811,PMIND,PMAXD,TMIND,TMAXD
DO 505 K=1,NPT
IF (PHI(K).GT.PI) GO TO 510
505 CONTINUE
GO TO 520
510 CONTINUE
DO 515 K=1,NPT
PHI(K)=PHI(K) - PI
515 THETA(K)= - THETA(K)
GO TO 570
520 CONTINUE
DO 525 K=1,NPT
IF (PHI(K).LE.(-PI)) GO TO 530
525 CONTINUE
GO TO 540
530 CONTINUE
DO 535 K=1,NPT
PHI(K)=PHI(K) + PI
535 THETA(K)= -THETA(K)
GO TO 570
540 CONTINUE
IPOS=0
INEG=0
DO 545 K=1,NPT
IF (PHI(K).GT.0.) IPOS=1
IF (PHI(K).LT.0.) INEG=1
545 CONTINUE
IF ((IPOS.EQ.1).AND.(INEG.EQ.1)) GO TO 570
C SUBCASE A: THE CAP WHICH WE HAVE DOES NOT CROSS THE LINE PHI=0.
C THEREFORE, THE BOUNDING CURVE CONSISTS OF TWO CLOSED CURVES: THE
C ONE WHICH WE HAVE AND ITS ANTI-PODE. ISIGN IS A SWITCH USED IN
C DETERMINING THE ANTI-PODAL CURVE.
IF ((IPOS.EQ.1).AND.(INEG.EQ.0)) ISIGN= -1
IF ((IPOS.EQ.0).AND.(INEG.EQ.1)) ISIGN= 1
IF ((IPOS.NE.1).AND.(INEG.NE.1)) THEN
JER=1
PRINT *,'ERROR IN BOUNDC (ICASE=5, SUBCASE A):'
PRINT *,'THE CURVE WE HAVE LIES ON THE LINE PHI=0, WHICH'
PRINT *,'IS VERY UNLIKELY.'
RETURN
END IF
IF (IND.EQ.1) THEN
WRITE (30,805) NPT
805 FORMAT (I4,' 2')
DO 560 K=1,NPT
BLONG=PHI(K)*180./PI
BLAT=THETA(K)*180./PI
560 WRITE (30,802) BLONG,BLAT
WRITE (30,805) NPT
DO 565 K=1,NPT
BLONG=PHI(K)*180./PI + ISIGN*180.
BLAT= -THETA(K)*180./PI
565 WRITE (30,802) BLONG,BLAT
CLOSE (30)
PRINT 815
815 FORMAT (' ',' NEITHER OF THE HOLES STRADDLES THE LINE:',
&' AXIS'/' ',' LONGITUDE = 180 DEGREES. HENCE, EACH HOLE IS',
&' BOUNDED BY A'/' ',' CLOSED CURVE.')
PRINT 715,NPT,BFILE
715 FORMAT (' ',' A SEQUENCE OF ',I3,' POINTS AROUND EACH OF',
&' THE CLOSED'/' ',' CURVES HAS BEEN WRITTEN TO THE FILE ',
&A10,'.')
END IF
RETURN
C SUBCASE B: THE CAP WHICH WE HAVE CROSSES THE LINE PHI=0. THEREFORE,
C THE ANTI-PODAL CAP IS BROKEN INTO TWO PIECES BY THE LINE PHI= +/-PI.
C THE BOUNDING CURVE CONSISTS OF THREE PIECES: ONE CLOSED CURVE AND TWO
C OPEN CURVES.
570 CONTINUE
IF (PHI(1).GE.0.) ISW= 1
IF (PHI(1).LT.0.) ISW= -1
DO 575 K=2,NPT
IF (PHI(K).GE.0.) KSW= 1
IF (PHI(K).LT.0.) KSW= -1
IF (KSW*ISW.LT.0) GO TO 580
575 CONTINUE
JER=1
PRINT *,'ERROR IN BOUNDC (ICASE=5, SUBCASE B):'
PRINT *,'IMPOSSIBLE EXIT FROM LOOP, STATEMENT 575.'
RETURN
580 K1=K-1
DO 585 K=K1+1,NPT
IF (PHI(K).GE.0.) KSW= 1
IF (PHI(K).LT.0.) KSW= -1
IF (KSW*ISW.GT.0) GO TO 590
585 CONTINUE
JER=1
PRINT *,'ERROR IN BOUNDC (ICASE=5, SUBCASE B):'
PRINT *,'IMPOSSIBLE EXIT FROM LOOP, STATEMENT 585.'
RETURN
590 K2=K-1
N1=K2 - K1
N2=NPT - 1 -N1
N3=NPT - 1 - K2
IF (IND.EQ.1) THEN
WRITE (30,805) NPT
DO 625 K=1,NPT
BLONG=PHI(K)*180./PI
BLAT=THETA(K)*180./PI
625 WRITE (30,802) BLONG,BLAT
WRITE (30,803) N1
DO 630 L=1,N1
K=K1 + L
BLONG=PHI(K)*180./PI + ISW*180.
BLAT= -THETA(K)*180./PI
630 WRITE (30,802) BLONG,BLAT
WRITE (30,803) N2
IF (N3.GT.0) THEN
DO 635 L=1,N3
K=K2 + L
BLONG=PHI(K)*180./PI - ISW*180.
BLAT= -THETA(K)*180./PI
635 WRITE (30,802) BLONG,BLAT
END IF
DO 645 L=N3+1,N2
K=L - N3
BLONG=PHI(K)*180./PI - ISW*180.
BLAT= -THETA(K)*180./PI
645 WRITE (30,802) BLONG,BLAT
CLOSE (30)
PRINT 816
PRINT 716,NPT,BFILE,N1,N2
END IF
816 FORMAT (' ',' ONE OF THE HOLES STRADDLES THE LINE: AXIS',
&' LONGITUDE ='/' ',' 180 DEGREES. HENCE, THIS HOLE IS BOUNDED',
&' BY TWO OPEN CURVES'/' ',' (ONE NEAR AXIS LONGITUDE = 180',
&' DEGREES, THE OTHER NEAR'/' ',' AXIS LONGITUDE = -180',
&' DEGREES). THE OTHER HOLE IS BOUNDED'/' ',' BY A CLOSED',
&' CURVE.')
716 FORMAT (' ',' A SEQUENCE OF ',I3,' POINTS AROUND THE CLOSED',
&' CURVE'/' ',' HAVE BEEN WRITTEN TO THE FILE ',A10,'; THESE ARE',
&' FOLLOWED'/' ',' BY ',I3,' POINTS ONE ONE OPEN CURVE AND ',I3,
&' POINTS ON THE'/' ',' OTHER.')
RETURN
END
C
C *** GRID ***
C
C SUBROUTINE TO CALCULATE THE MIN. AND MAX. VALUES OF RHO (ANGLE OF
C ROTATION) ON A RECTANGULAR GRID IN THE AXIS LONGITUDE-AXIS LATITUDE
C PLANE.
C INPUTS:
C ALAT1 - LATITUDE VALUE FOR FIRST ROW
C ALAT2 - LATITUDE VALUE FOR LAST ROW
C ALONG1 - LONGITUDE VALUE FOR FIRST ROW
C ALONG2 - LONGITUDE VALUE FOR LAST ROW
C (THESE VALUES ARE IN DEGREES.)
C QF - THE MATRIX QTILDE IN FULL STORAGE MODE
C MF - THE MATRIX M IN FULL STORAGE MODE
C W - ARRAY CONTAINING THE EIGENVECTORS OF M
C ICASE - INDICATOR FOR THE TYPE OF SET A
C IND - INDICATOR
C IND=0: BOX OUT THE CONFIDENCE REGION
C IND=1: WRITE UPPER AND LOWER SURFACES TO FILES 31 AND 32, RESP.
C OUTPUTS:
C JER - ERROR INDICATOR
C 0 - ALL IS WELL
C 1 - AN ERROR HAS OCCURRED
C DEFINED BY PARAMETER STATEMENTS:
C NLAT - NO. OF ROWS IN THE GRID
C NLONG - NO. OF COLS. IN THE GRID
C IF IND=1, THE SUBROUTINE GRID WRITES THE RHO VALUES FOR THE UPPER
C SURFACE TO FILE 31; RHO VALUES FOR THE LOWER SURFACE TO FILE 32.
C (AN EXCEPTION OCCURS WHEN ICASE=7: IN THIS CASE ONLY FILE 31 IS
C USED; THE MIN. VALUE OF RHO IS ZERO FOR EACH AXIS.)
SUBROUTINE GRID(ALAT1,ALAT2,ALONG1,ALONG2,QF,MF,W,ICASE,IND,JER)
PARAMETER (NLAT=101)
PARAMETER (NLONG=201)
REAL QF(4,4),MF(3,3),W(3,3)
REAL U(3),MLIM,MTU(3),AF(2,2),BU(3)
REAL MU(2),V(2,2),WORK(4)
REAL DRHO(NLONG),CRHO(NLONG)
CHARACTER UFILE*10,LFILE*10
DATA PI/3.1415927/
JER=0
IF (IND.EQ.1) THEN
PRINT *,' TYPE NAME OF THE UPPER SURFACE FILE, I.E., THE FILE'
PRINT *,' TO RECEIVE MAX. ANGLE OF ROTATION VALUES. (MAX.='
PRINT *,' 10 CHARACTERS.)'
READ 999,UFILE
999 FORMAT (A)
OPEN (31,FILE=UFILE,STATUS='UNKNOWN')
REWIND 31
IF (ICASE.LT.7) THEN
PRINT *,' TYPE NAME OF THE LOWER SURFACE FILE, FOR MIN.'
PRINT *,' ANGLE OF ROTATION VALUES. (MAX.=10 CHARACTERS.)'
READ 999,LFILE
OPEN (32,FILE=LFILE,STATUS='UNKNOWN')
REWIND 32
END IF
END IF
DLONG=(ALONG2 - ALONG1)/(NLONG - 1.0)
DLAT=(ALAT2 - ALAT1)/(NLAT - 1.0)
C IRHO IS A SWITCH FOR CALCULATION OF MIN. AND MAX. VALUES OF RHO.
C IRHO= 0: NO VALUES OF RHO HAVE BEEN CALCULATED
C IRHO= 1: AT LEAST ONE VALUE OF RHO HAS BEEN CALCULATED
C WHEN THE FIRST VALUE OF RHO IS CALCULATED, IRHO IS SWITCHED FROM 0
C TO 1.
IRHO=0
C OUTER LOOP: ON THE LATITUDE VALUES. THE ARRAYS CRHO AND DRHO
C CONTAIN THE VALUES OF RHO FOR THE LOWER AND UPPER SURFACES, RESP.,
C ALONG THE CURRENT ROW. THESE ARRAYS ARE INITIALLY SET TO THE
C VALUE -1000000.; THEN THE CORRECT VALUES FOR GRID POINTS WHICH
C COORESPOND TO ADMISSIBLE AXES ARE FILLED IN. THE VALUES
C -1000000. WHICH ARE LEFT SIGNAL NON-ADMISSIBLE AXES.
DO 400 I=1,NLAT
CLAT=ALAT1 + (I-1)*DLAT
DO 50 J=1,NLONG
DRHO(J)= -1000000.
50 CRHO(J)= -1000000.
C INNER LOOP: ON THE LONGITUDE VALUES. IF ICASE.LE.5, AN AXIS U IS
C SKIPPED IF TRANSPOSE(U)*M*U.GT.(AZERO - 1). IN ADDITION, IF ICASE
C .LE.3, U IS SKIPPED IF IT IS NOT IN THE DESIRED CAP. FOR A PARTIC-
C ULAR AXIS U WHICH IS ACCEPTED, THE PERMISSIBLE RANGE OF ANGLES OF
C ROTATION IS CALCULATED BASED ON THE MATRIX A(U).
DO 350 J=1,NLONG
CLONG=ALONG1 + (J-1)*DLONG
CPHI=CLONG*PI/180.
CTHETA=CLAT*PI/180.
COST=COS(CTHETA)
U(1)=COST*COS(CPHI)
U(2)=COST*SIN(CPHI)
U(3)=SIN(CTHETA)
IF (ICASE.LE.5) THEN
AZERO=QF(1,1)
MLIM=AZERO - 1.0
DO 155 K=1,3
MTU(K)=0.0
DO 155 L=1,3
155 MTU(K)=MTU(K) + MF(K,L)*U(L)
UMU=0.0
DO 165 K=1,3
165 UMU=UMU + MTU(K)*U(K)
IF (UMU.GT.MLIM) GO TO 350
IF (ICASE.LE.3) THEN
DOTP=0.0
DO 190 K=1,3
190 DOTP=DOTP + U(K)*W(K,3)
IF (DOTP.LE.0.0) GO TO 350
END IF
END IF
AF(1,1)=QF(1,1)
SUM=0.0
DO 270 K=1,3
270 SUM=SUM + QF(K+1,1)*U(K)
AF(1,2)=SUM
AF(2,1)=AF(1,2)
DO 275 K=1,3
BU(K)=0.0
DO 275 L=1,3
275 BU(K)=BU(K) + QF(K+1,L+1)*U(L)
SUM=0.0
DO 285 K=1,3
285 SUM=SUM + BU(K)*U(K)
AF(2,2)=SUM
CALL JACOB2(AF,2,2,MU,V,2,WORK,NROT)
IF (NROT.LE.0) THEN
JER=1
PRINT *,'ERROR IN GRID, ON A CALL TO JACOB2.'
PRINT *,'ROW = ',I,' COLUMN = ',J
PRINT *,'NROT = ', NROT
RETURN
END IF
C IF ICASE.LT.7, ARRANGE THAT V(2,1).GE.0., SO THAT RHOSTAR COMES OUT IN
C THE RANGE FROM 0 TO 360 DEGREES.
IF (ICASE.LT.7) THEN
IF (V(2,1).LT.0.) THEN
DO 302 K=1,2
302 V(K,1)= -V(K,1)
END IF
END IF
C IF ICASE.EQ.7, ARRANGE THAT V(1,1).GE.0., SO THAT RHOSTAR IS IN THE
C RANGE FROM -180 TO 180 DEGREES.
IF (ICASE.EQ.7) THEN
IF (V(1,1).LT.0.) THEN
DO 308 K=1,2
308 V(K,1)= -V(K,1)
END IF
END IF
ANGLE=ATAN2(V(2,1),V(1,1))
RHOSTAR=(2.*ANGLE)*180./PI
IF (MU(1).GE.1.0) THEN
RHOINC=0.0
ELSE IF (MU(2).LE.1.0) THEN
RHOINC=180.
ELSE
TD=SQRT((1.0 - MU(1))/(MU(2) - 1.0))
DELTA=ATAN(TD)
RHOINC=(2.*DELTA)*180./PI
END IF
DRHO(J)=RHOSTAR + RHOINC
CRHO(J)=RHOSTAR - RHOINC
IF (IRHO.NE.1) THEN
IRHO=1
RMIND=CRHO(J)
RMAXD=DRHO(J)
ELSE
IF (CRHO(J).LT.RMIND) RMIND=CRHO(J)
IF (DRHO(J).GT.RMAXD) RMAXD=DRHO(J)
END IF
350 CONTINUE
C END OF INNER LOOP.
IF (IND.EQ.1) THEN
WRITE (31,905) (DRHO(J),J=1,NLONG)
IF (ICASE.LT.7) WRITE (32,905) (CRHO(J),J=1,NLONG)
905 FORMAT (4E16.8)
END IF
400 CONTINUE
C END OF OUTER LOOP.
IF (IRHO.EQ.0) THEN
PRINT *,'UNABLE TO CALCULATE MIN. AND MAX. VALUES OF THE ANGLE'
PRINT *,'OF ROTATION, SINCE NONE OF THE GRID POINTS CORRESPOND'
PRINT *,'TO ADMISSIBLE AXES OF ROTATION.'
PRINT 951,ALONG1,ALONG2,ALAT1,ALAT2,NLONG,NLAT
RETURN
END IF
951 FORMAT (' ','MIN., MAX. AXIS LONGITUDE OVER GRID: ',F7.1,',',
&F7.1,' DEGREES.'/' ','MIN., MAX. AXIS LATITUDE OVER GRID: ',
&F7.1,',',F7.1,' DEGREES.'/' ','NO. OF LONGITUDE VALUES, I.E.,',
&' COLS., IN GRID= ',I3/' ','NO. OF LATITUDE VALUES, I.E.,',
&' ROWS, IN GRID= ',I3)
IF (ICASE.EQ.7) RMIND=0.0
PRINT 952,RMIND,RMAXD
952 FORMAT (' ','MIN., MAX. ANGLE OF ROTATION OVER THE CONFIDENCE',
&' REGION:'/' ',5X,F10.5,',',F10.5,' DEGREES.'/' ','NOTE: THESE',
&' VALUES ARE THE MIN. AND MAX. OVER A RECTANG-'/' ','ULAR GRID',
&' SUPERIMPOSED ON THE SET OF ADMISSIBLE AXES.')
PRINT 951,ALONG1,ALONG2,ALAT1,ALAT2,NLONG,NLAT
IF (IND.EQ.1) THEN
IF (ICASE.LT.7) THEN
CLOSE (31)
CLOSE (32)
PRINT 916,UFILE,NLONG,LFILE
ELSE IF (ICASE.EQ.7) THEN
CLOSE (31)
PRINT 917,UFILE,NLONG
END IF
END IF
916 FORMAT (' ',' MAX. ANGLE OF ROTATION VALUES HAVE BEEN',
&' WRITTEN TO'/' ',' THE FILE ',A10,' (BY GRID ROWS, WITH ',I3,
&' VALUES'/' ',' TO A ROW). MIN. ANGLE OF ROTATION VALUES',
&' (ALSO BY ROW)'/' ',' HAVE BEEN WRITTEN TO THE FILE ',A10,'.')
917 FORMAT (' ',' MAX. ANGLE OF ROTATION VALUES HAVE BEEN',
&' WRITTEN TO'/' ',' THE FILE ',A10,' (BY GRID ROWS, WITH ',I3,
&' VALUES'/' ',' TO A ROW). NOTE: ALL MIN. ANGLE OF ROTATION',
&' VALUES ARE'/' ','ZERO.')
RETURN
END
C
C *** MEIG ***
C
C SUBROUTINE TO CALCULATE THE MATRIX:
C M = (AZERO - 1)*B - DVECT*TRANSPOSE(DVECT)
C AND TO FIND ITS EIGENVALUES AND EIGENVECTORS.
C INPUTS:
C QF - THE MATRIX QTILDE IN FULL STORAGE MODE
C UXX - OPTIMAL AXIS OF ROTATION
C OUTPUTS:
C NU - A 3-VECTOR CONTAINING THE EIGENVALUES OF M
C W - A 3X3 MATRIX CONTAINING (AS ITS COLUMNS) THE EIGENVECTORS OF M
C MF - THE MATRIX M (IN FULL STORAGE MODE)
C ICASE - INDICATOR FOR THE TYPE OF REGION
C JER - AN ERROR INDICATOR
C 0 - ALL IS WELL
C 1 - AN ERROR HAS OCCURRED
SUBROUTINE MEIG(QF,UXX,NU,W,MF,ICASE,JER)
REAL QF(4,4),MF(3,3),DVECT(3),B(3,3)
REAL NU(3),W(3,3),D(3),Z(3,3),WORK(6)
REAL UXX(3)
INTEGER KPOS(3)
JER=0
AZERO=QF(1,1)
DO 15 I=1,3
15 DVECT(I)=QF(I+1,1)
DO 18 I=1,3
DO 18 J=1,3
18 B(I,J)=QF(I+1,J+1)
DO 30 I=1,3
DO 30 J=1,3
30 MF(I,J)=(AZERO - 1.0)*B(I,J) - DVECT(I)*DVECT(J)
CALL JACOB2(MF,3,3,D,Z,3,WORK,NROT)
IF (NROT.LE.0) THEN
JER=1
PRINT *,'ERROR IN MEIG, ON A CALL TO JACOB2:'
PRINT *,'NROT = ',NROT
RETURN
END IF
C CHECK THAT M HAS ONE NEGATIVE AND TWO POSITIVE EIGENVALUES.
NNEG=0
NPOS=0
DO 40 I=1,3
IF (D(I).LT.0.) THEN
NNEG=NNEG + 1
KNEG=I
ELSE IF (D(I).GT.0.) THEN
NPOS=NPOS + 1
KPOS(NPOS)=I
END IF
40 CONTINUE
IF ((NNEG.NE.1).OR.(NPOS.NE.2)) THEN
JER=1
PRINT *,'ERROR IN MEIG:'
PRINT *,'IT IS NOT THE CASE THAT M HAS ONE NEGATIVE EIGENVALUE'
PRINT *,'AND TWO POSITIVE ONES. THERE IS SOME ERROR, AS M'
PRINT *,'SHOULD HAVE EIGENVALUES OF THIS TYPE.'
RETURN
END IF
C COUNT THE NO. OF POSITIVE EIGENVALUES OF M WHICH ARE GREATER THAN
C AZERO - 1.
NLARG=0
DO 47 J=1,2
K=KPOS(J)
IF (D(K).GT.(AZERO - 1.0)) NLARG=NLARG + 1
47 CONTINUE
IF (NLARG.EQ.2) GO TO 100
IF (NLARG.EQ.1) GO TO 200
C NLARG=0. ANY AXIS IS ADMISSIBLE.
ICASE=6
RETURN
C NLARG=2. THE SET OF ADMISSIBLE AXES CONSISTS OF TWO (ANTI-PODAL)
C CAPS. SET NU(3)=THE NEGATIVE EIGENVALUE AND ARRANGE THAT THE EIGEN-
C VECTOR CORRESP. TO NU(3) MAKES AN ACUTE ANGLE WITH THE OPTIMAL AXIS.
C IF THE PREFERRED CAP DOES NOT CONTAIN ONE OF THE POLES, ICASE=1.
C OTHERWISE, ICASE=2 OR 3, DEPENDING ON WHETHER THE POLE INVOLVED IS
C THE NORTH OR THE SOUTH.
100 CONTINUE
NU(3)=D(KNEG)
DO 147 I=1,3
147 W(I,3)=Z(I,KNEG)
DO 150 J=1,2
K=KPOS(J)
NU(J)=D(K)
DO 150 I=1,3
150 W(I,J)=Z(I,K)
DOTP=0.0
DO 157 I=1,3
157 DOTP=DOTP + W(I,3)*UXX(I)
IF (DOTP.LT.0.) THEN
DO 160 I=1,3
160 W(I,3)= -W(I,3)
END IF
IF (MF(3,3).GT.(AZERO - 1.0)) THEN
ICASE=1
GO TO 300
ELSE
IF (W(3,3).GE.0.) ICASE=2
IF (W(3,3).LT.0.) ICASE=3
GO TO 300
END IF
C NLARG=1. THE SET OF ADMISSIBLE AXES IS THE COMPLEMENT OF TWO
C EXCLUDED (ANTI-PODAL) CAPS. SET NU(1)=THE NEGATIVE EIGENVALUE.
C NOTE: NU(2).LE.NU(3), SO NU(3) IS THE POSITIVE EIGENVALUE WHICH
C IS GREATER THAN AZERO - 1. IF THE EXCLUDED CAPS CONTAIN THE
C POLES, ICASE=4. OTHERWISE, ICASE=5. IF ICASE=4, ARRANGE THAT
C THE EIGENVECTOR CORRESP. TO NU(3) IS IN THE SAME CAP AS THE
C NORTH POLE.
200 CONTINUE
NU(1)=D(KNEG)
DO 205 I=1,3
205 W(I,1)=Z(I,KNEG)
DO 210 J=2,3
K=KPOS(J-1)
NU(J)=D(K)
DO 210 I=1,3
210 W(I,J)=Z(I,K)
IF (MF(3,3).LE.(AZERO - 1.0)) THEN
ICASE=5
GO TO 300
ELSE
ICASE=4
IF (W(3,3).LT.0.) THEN
DO 235 I=1,3
235 W(I,3)= -W(I,3)
END IF
END IF
C MAKE SURE THAT THE COLUMNS OF W FORM A RIGHT-HANDED SYSTEM.
300 CONTINUE
DET= W(1,1)*W(2,2)*W(3,3)
& + W(1,2)*W(2,3)*W(3,1)
& + W(1,3)*W(3,2)*W(2,1)
& - W(1,3)*W(2,2)*W(3,1)
& - W(1,2)*W(2,1)*W(3,3)
& - W(1,1)*W(2,3)*W(3,2)
IF (DET.LE.0.) THEN
DO 305 I=1,3
305 W(I,1)= -W(I,1)
END IF
RETURN
END
C
C *** EVALF ***
C
C SUBROUTINE TO EVALUATE THE FUNCTION
C F(PHIT) = 1./SQRT((DPHI/DPHIT)**2 + (DTHETA/DPHIT)**2)
C NOTE: F(PHIT) = DPHIT/DS, WHERE S IS THE ARC LENGTH ALONG THE
C BOUNDING CURVE (MEASURED IN RADIANS).
C ADDITIONAL NOTE: IN THE PROCESS OF EVALUATING F(PHIT), EVALF
C ALSO CALCULATES PHI, THETA, AND THE VECTOR U.
C INPUTS:
C PHIT - PARAMETER ALONG THE BOUNDING CURVE
C NU - ARRAY CONTAINING THE EIGENVALUES OF THE MATRIX M
C W - ARRAY CONTAINING THE EIGENVECTORS OF THE MATRIX M
C AZERO - QF(1,1)
C ICODE - INDICATOR FOR MODE OF USAGE
C 0 - FIRST POINT ON BOUNDING CURVE
C POS. - SUBSEQUENT POINTS
C OLDU - THE U -VECTOR FOR THE IMMEDIATELY PRECEDING POINT, IN THE
C CASE ICODE IS POSITIVE
C OPHI - THE PHI VALUE FOR THE IMMEDIATELY PRECEDING POINT, IN THE
C CASE ICODE IS POSITIVE
C OUTPUTS:
C FVAL - VALUE OF THE FUNCTION F(PHIT)
C PHI - LONGITUDE CORRESP. TO PHIT
C THETA - LATITUDE CORRESP. TO PHIT
C U - AXIS OF ROTATION CORRESP. TO PHIT
C JER - ERROR INDICATOR
C 0 - ALL IS WELL
C 1 - AN ERROR HAS OCCURRED
C PHIT AND THETAT ARE THE LONGITUDE AND LATITUDE, RESPECTIVELY, MEASURED
C IN THE W-COORDINATE SYSTEM (I.E., THE COORDINATE SYSTEM DETERMINED BY
C THE EIGENVECTORS OF M).
SUBROUTINE EVALF(PHIT,NU,W,AZERO,FVAL,PHI,THETA,U,ICODE,
&OLDU,OPHI,JER)
REAL NU(3),W(3,3),U(3),ETA(3),DUDPT(3)
REAL DEDPT(3),LAMBDA
REAL NUM,DENOM,OLDU(3)
DATA PI/3.1415927/
JER=0
COSPT=COS(PHIT)
SINPT=SIN(PHIT)
DENOM=NU(1)*COSPT*COSPT + NU(2)*SINPT*SINPT - NU(3)
NUM=AZERO - 1.0 - NU(3)
COSTT=SQRT(NUM/DENOM)
THETAT=ACOS(COSTT)
ETA(1)=COSTT*COSPT
ETA(2)=COSTT*SINPT
ETA(3)=SIN(THETAT)
DO 10 K=1,3
U(K)=0.0
DO 10 L=1,3
10 U(K)=U(K) + W(K,L)*ETA(L)
LAMBDA=(NU(2) - NU(1))/NUM
DEDPT(1)= -ETA(2)*(LAMBDA*ETA(1)*ETA(1) + 1.)
DEDPT(2)= -ETA(1)*(LAMBDA*ETA(2)*ETA(2) - 1.)
SS=ETA(1)*ETA(1) + ETA(2)*ETA(2)
DEDPT(3)=LAMBDA*ETA(1)*ETA(2)*SS/ETA(3)
DO 20 K=1,3
DUDPT(K)=0.0
DO 20 L=1,3
20 DUDPT(K)=DUDPT(K) + W(K,L)*DEDPT(L)
COST2=U(1)*U(1) + U(2)*U(2)
DPDPT=( -U(2)*DUDPT(1) + U(1)*DUDPT(2))/COST2
DTDPT=DUDPT(3)/SQRT(COST2)
FVAL=1./SQRT(DPDPT*DPDPT + DTDPT*DTDPT)
PHI=ATAN2(U(2),U(1))
THETA=ASIN(U(3))
C MODIFY PHI IF NECESSARY.
IF (ICODE.EQ.0) THEN
RETURN
END IF
ALPHA=U(1)*OLDU(1) + U(2)*OLDU(2)
BETA= -U(1)*OLDU(2) + U(2)*OLDU(1)
DELPHI=ATAN2(BETA,ALPHA)
EPHI=OPHI + DELPHI
IF (ABS(EPHI - PHI).LT.PI/180.) THEN
RETURN
ELSE IF (ABS(EPHI - PHI - 2.*PI).LT.PI/180.) THEN
PHI=PHI + 2.*PI
RETURN
ELSE IF (ABS(EPHI - PHI + 2.*PI).LT.PI/180.) THEN
PHI=PHI - 2.*PI
RETURN
ELSE
JER=1
PRINT *,'ERROR IN EVALF: UNABLE TO DETERMINE PHI, FOR A'
PRINT *,'POINT ON THE BOUNDING CURVE.'
RETURN
END IF
END
C
C *** LARGIN ***
C
C LARGIN(X) = THE LARGEST INTEGER WHICH IS LESS THAN OR EQUAL TO X.
INTEGER FUNCTION LARGIN(X)
LARGIN=INT(X)
IF (X.LT.0.) THEN
R=X - FLOAT(LARGIN)
IF (R.LT.0.) LARGIN=LARGIN -1
END IF
RETURN
END