C C THE PROGRAM (KRNLS) COMPUTES AND PLOTS A SERIES OF C C S C K AND K KERNLES (N,M=0...NMAX) FOR A GIVEN DIRECTION "AZ". C NM NM C WHERE C C C K =COS(M*AZ)*F (PSI0,PSI), C NM NM C S C K =SIN(M*AZ)*F (PSI0,PSI). WHERE C NM NM C C F =K*INTERAL OF (F1*F2*F3) C NM C WHERE C _ INF _ _ C F1=SIN(K*PSI) ,F2=P (K*PSIS) ,F3=ZIG(P (PSI)*P (PSIS)/(I-1)/C). C NM I=M IM IM C WHERE C=2 FOR M=0 AND C=4 FOR M>0. C C C INPUT DATA- C NMAX- IS THE MAXIMUM VALUE FOR N. C NCAP- MAXIMUM DEGREES UP TO WHICH THE ASSOCIATED C LEGENDRE FUNCTIONS ARE COMPUTED. C PSI0- SPHERICAL RADIUS OF THE CAP(AREA),IN DEGREES. C INFINT- THE MAXIMUM LIMIT OF THE SUMMATION IN F3. C MPOINT- THE NUMBER OF BASE POINTS USED BY THE C GAUSS LEGENDRE QUADRATURE.IT CAN TAKE ON 2,3,4,5,6 C 10 OR 15. C DSDIS- AN INCREMENT OF SPHERICAL DISTANCE,IN DEGREES. C AZ- AN AZIMUTH (DIRECTION),IN DEGREES. C NOTE: THE INPUT VARIABLES ARE ASSIGNED IN THE C PROGRAM. C OUTPUT DATA- C THE GRAPH OF THE KERNELS. C M. NAJAFI 9/1980 IMPLICIT REAL*8(A-H,O-Z) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CZ/ ZPP(20,15),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M INTEGER C(81) REAL*8 PN(16,16),KRNL(60,81),X(60),A(81),FNM(20,20) EXTERNAL F IPN=16 IPNM=541 PI=DARCOS(-1.D0) DR=PI/180.D0 C NCAP=10 C INFINT=180 C MPOINT=15 C C NMAX CAN TAKE ON VALUES FROM 0 UP TO NCAP-1. C NMAX=5 C PSI0=10.D0 C DSDIS=1.D0 C AZ=0.D0 C NMP1=NMAX+1 C CONSTRUCTING A TABLE OF NORMALIZED ASSOCIATED LEGENDRE FUNCTIONS C UP TO DEGREE AND ORDER NCAP FOR ANGELS FROM 0 TO 90 DEGREES IN C STEPS OF 10 MINUTES. C NCAP1=NCAP+1 SDIS=0.D0 PSINT=10.D0/60 DO 10 IR=1,IPNM PHI=90.D0-SDIS CALL LEGDRE(1,PHI,PN,IPN,NCAP1) K=0 DO 11 I=1,NCAP1 DO 11 J=1,I K=K+1 11 PNM(IR,K)=PN(I,J) 10 SDIS=SDIS+PSINT JPNM=K C C COMPUTING THE KERNELS FOR A DIRECTION (AZ=...) DEGREES. C AZ=AZ*DR I=1 C(I)=100 A(I)=1.D0 IF(NMAX.EQ.0) GO TO 100 DO 20 N=1,NMAX I=I+1 C(I)=100+N*10 A(I)=1.D0 DO 21 M=1,N I=I+1 C(I)=100+N*10+M A(I)=DCOS(M*AZ) I=I+1 C(I)=200+N*10+M A(I)=DSIN(M*AZ) 21 CONTINUE 20 CONTINUE 100 CONTINUE PSI=0.D0 IP=0 30 CONTINUE IP=IP+1 X(IP)=PSI CALL ZGMAPP(MPOINT,NMAX,0.D0,PSI0*DR) AA=0.D0 BB=PSI0*DR CALL ISOKNL(AA,BB,NMAX,MPOINT,FNM,II,JJ) JP=1 KRNL(IP,JP)=A(JP)*FNM(1,1) IF(NMAX.EQ.0) GO TO 200 DO 31 N=1,NMAX JP=JP+1 KRNL(IP,JP)=A(JP)*FNM(N+1,1) DO 32 M=1,N JP=JP+1 KRNL(IP,JP)=A(JP)*FNM(N+1,M+1) JP=JP+1 KRNL(IP,JP)=A(JP)*FNM(N+1,M+1) 32 CONTINUE 31 CONTINUE 200 CONTINUE IF(PSI.GT.PSI0) DSDIS=10.D0 PSI=PSI+DSDIS IF(PSI.LE.180.D0) GO TO 30 C C PLOTTING THE KERNELS. C DO 40 J=1,JP DO 41 I=1,IP 41 A(I)=KRNL(I,J) CALL GRAPH(X,A,IP,' SPH.DIS',' KERNEL') PRINT,' ' PRINT,' ' PRINT,' & (',J,')' IC=C(J)/10.D0 IF(C(J).GT.200) GO TO 18 I1=IC-10 I2=C(J)-10*IC WRITE(6,14) I1,I2 GO TO 40 18 I1=IC-20 I2=C(J)-10*IC WRITE(6,15) I1,I2 40 CONTINUE 14 FORMAT('0',83X,'A',I1,I1) 15 FORMAT('0',83X,'B',I1,I1) RETURN END SUBROUTINE ISOKNL(A,B,INMAX,MPOINT,FNM,IROW,ICUL) C C ISOKNL COMPUTES THE F (PSI0,PSI) FUNCTIONS OF (SEE MAIN PROGRAM) C NM C OF N,M=0...INMAX FOR A GIVEN SPHERICAL DISTANCE PSI. C INPUTE DATA: A,B- LOWER AND UPPER LIMIT OF THE INTEGRATION C IN RADIANS. C INMAX- THE MAXIMUM LIMIT OF N. C MPOINT- NUMBER OF BASE POINTS USED BY ROUTINE GAUSS. C OUTPUT DATA: FNM(IROW,ICUL)- THE ARRAY CONTAINING THE FUNCTION C VALUES. C M. NAJAFI, 10/1980 IMPLICIT REAL*8 (A-H,O-Z) REAL*8 FNM(20,20) EXTERNAL F COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/PP/ PSI0,PSI,N,M COMMON/CZ/ ZPP(20,15),INDEX,IZP,JZP,INFINT SCALE=180.D0/PSI0 PI=DARCOS(-1.D0) DR=PI/180.D0 DO 150 I=1,20 DO 150 J=1,20 150 FNM(I,J)=0.D0 C.....PSI AND PSI0 HAVE TO BE GIVEN IN DEGREES. NBP1=INMAX+1 DO 100 IROW=1,NBP1 N=IROW-1 DO 200 ICUL=1,IROW M=ICUL-1 FNM(IROW,ICUL)=GAUSS(A,B,MPOINT,F)*SCALE 200 CONTINUE 100 CONTINUE RETURN END SUBROUTINE ZGMAPP(MP,MMAX,A,B) C C ZGMAPP COMPUTES THE F3 FUNCTION (SEE MAIN PROGRAM) OF M=0...MMAX C FOR A GIVEN PSI AND FOR THE SAME BASE POINTS USED BY THE C ROUTINE GAUSS. C INPUT DATA: MP- NUMBER OF THE BASE POINTS WHICH INTRODUCES C THE BASE POINTS TO THE ROUTINE. C MMAX- THE MAXIMUM ORDER (M). C A,B- LOWER AND UPPER LIMIT OF THE INTEGRATION. C OUTPUT DATA: ZPP(IZP,JZP)- AN ARRAY CONTAINING THE VALUES C OF F3 FUNCTION. C M. NAJAFI 10/1980 IMPLICIT REAL*8 (A-H,O-Z) REAL*8 ZI(32),XI(41) INTEGER NPOINT(8),KEY(9) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CZ/ ZPP(20,15),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M C C.....PRESET NPOINT,KEY,ZI,AND WEIGHT ARRAYS..... C DATA NPOINT/2,3,4,5,6,10,15,16/ C DATA KEY/1,2,4,6,9,12,17,25,33/ C DATA(ZI(I),I=1,24)/0.577350269D0,0.0D0,0.774596669D0,0.339981044 &D0, 0.861136312D0,0.0D0 ,0.538469310D0,0.906179846D0, & 0.238619186D0,0.661209387D0,0.932469514D0,0.148874339D0, & 0.433395394D0,0.679409568D0,0.865063367D0,0.973906529D0, &0.0D0,0.201194093997D0,0.394151347078D0,0.570972172609D0,0.72441 &7731360D0,0.848206583410D0,0.937273392401D0,0.987992518020D0/ DATA (ZI(I),I=25,32)/0.09501250979509,0.28160355077652,0.458016777 &67026,0.61787624436991,0.75540440831005,0.86563120242217,0.9445750 &2305806,0.98940093496853/ PI=DARCOS(-1.D0) DR=PI/180.D0 C C.....FIND SUBSCRIPT OF FIRST ZI AND WEIGHT VALUE...... C DO 1 I=1,8 IF(MP.EQ.NPOINT(I)) GO TO 2 1 CONTINUE C C.....INVALID M USED. C WRITE(6,11) 11 FORMAT(' ','******INVALID M HAS BEEN USED.******') GAUSS=0.D0 RETURN C C.....SET UP INITIAL PARAMETERS...... C 2 JFIRST=KEY(I) JLAST=KEY(I+1)-1 C=(B-A)/2.D0 D=(B+A)/2.D0 INFT=INFINT JZP=MMAX+1 K=0 DO 5 J=JFIRST,JLAST K=K+1 IF(ZI(J).NE.0.D0) GO TO 6 XI(K)=D GO TO 5 6 CONTINUE XI(K)=D+ZI(J)*C K=K+1 XI(K)=D-ZI(J)*C 5 CONTINUE IZP=K DO 7 II=1,IZP PSISTR=XI(II) DO 8 JJ=1,JZP M=JJ-1 FOUR=4.D0 IF(M.EQ.0) FOUR=2.D0 I1=M IF(M.LT.2) I1=2 Z=0.D0 P12=0.D0 P22=0.D0 DO 10 I=I1,NCAP P11=P12 P21=P22 P12=P(I,M,PSI) P22=P(I,M,PSISTR/DR) CALL TRAPS(0,0,1000,0,0) Z1=P12*P22/(I-1)/FOUR 10 Z=Z+Z1 KKK=NCAP+1 C'''''''''''''''''''''''''''''' C P11=P(I1,M,PSI) C P21=P(I1,M,PSISTR/DR) C P12=P(I1+1,M,PSI) C P22=P(I1+1,M,PSISTR/DR) C CALL TRAPS(0,0,1000,0,0) C Z1=P11*P21/(I1-1)/FOUR C CALL TRAPS(0,0,1000,0,0) C Z2=P12*P22/I1/FOUR C Z=Z1+Z2 C KKK=I1+2 C'''''''''''''''''''''''''''''' COSDD=DCOS(PSI*DR) COSTR=DCOS(PSISTR) DO 30 I=KKK,INFT CO1=(2*I-1)*1.D0/(I-M)*H(I,M)/H(I-1,M) CO2=(I+M-1)*1.D0/(2*I-1)*H(I-1,M)/H(I-2,M) HOLD=P12 P12=CO1*(COSDD*P12-CO2*P11) P11=HOLD HOLD=P22 P22=CO1*(COSTR*P22-CO2*P21) P21=HOLD Z1=P12*P22/(I-1)/FOUR 30 Z=Z+Z1 ZPP(II,JJ)=Z 8 CONTINUE 7 CONTINUE RETURN END DOUBLE PRECISION FUNCTION GAUSS(A,B,MP,F) C C GAUSS-LEGENDRE QUADRATURE EVALUATES AN INTEGRAL BETWEEN THE C INTERVAL (A,B),SEE CARNAHAN 1976. C INPUT DATA: A,B- LOWER AND UPPER LIMIT OF THE INTEGRATION. C MP- DEGREE OF LEGENDRE POLYNOMIAL WHOSE ZEROES C ARE USED AS THE BASE POINTS IN THE ROUTINE. C F- VALUE OF THE INTEGRAND AT A POINT. C OUTPUT DATA: GAUSS- VALUE OF THE INTEGRAL. C M. NAJAFI, 10/1980 IMPLICIT REAL*8(A-H,O-Z) REAL*8 ZI(32),WEIGHT(32) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CZ/ ZPP(20,15),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M INTEGER NPOINT(8),KEY(9) C C.....PRESET NPOINT,KEY,ZI,AND WEIGHT ARRAYS..... C DATA NPOINT/2,3,4,5,6,10,15,16/ C DATA KEY/1,2,4,6,9,12,17,25,33/ C DATA(ZI(I),I=1,24)/0.577350269D0,0.D0,0.774596669D0,0.339981044 &D0, 0.861136312D0,0.0D0 ,0.538469310D0,0.906179846D0, & 0.238619186D0,0.661209387D0,0.932469514D0,0.148874339D0, & 0.433395394D0,0.679409568D0,0.865063367D0,0.973906529D0, &0.0D0,0.201194093997D0,0.394151347078D0,0.570972172609D0,0.72441 &7731360D0,0.848206583410D0,0.937273392401D0,0.987992518020D0/ DATA (ZI(I),I=25,32)/0.09501250979509,0.28160355077652,0.458016777 &67026,0.61787624436991,0.75540440831005,0.86563120242217,0.9445750 &2305806,0.98940093496853/ C DATA (WEIGHT(I),I=1,24)/1.D0,0.888888889D0,0.555555556D0,0.652145 &155D0, 0.347854845D0,0.568888889D0,0.478628671D0,0.236926885D0, & 0.467913935D0,0.360761573D0,0.171324492D0,0.295524225D0, &0.269266719D0,0.219086363D0,0.149451349D0,0.066671344D0,0.2025782 &41926D0,0.198431485327D0,0.186161000116D0,0.166269205817D0,0.1395 &70677926D0,0.107159220467D0,0.070366047488D0,0.030753241996D0/ DATA (WEIGHT(I),I=25,32)/0.18945061039604,0.18260341508958,0.16915 &651936061,0.14959598877491,0.12462897132867,0.09515851170809,0.062 &25352386958,0.02715245945325/ C C.....FIND SUBSCRIPT OF FIRST ZI AND WEIGHT VALUE...... C DO 1 I=1,8 IF(MP.EQ.NPOINT(I)) GO TO 2 1 CONTINUE C C.....INVALID M USED. C WRITE(6,10) 10 FORMAT(' ','******INVALID M HAS BEEN USED.******') GAUSS=0.D0 RETURN C C.....SET UP INITIAL PARAMETERS...... C 2 JFIRST=KEY(I) JLAST=KEY(I+1)-1 C=(B-A)/2.D0 D=(B+A)/2.D0 C C.....ACCUMULATE THE SUM IN THE M-POINT FORMULA..... SUM=0.D0 INDEX=0 DO 5 J=JFIRST,JLAST IF(ZI(J).NE.0.D0) GO TO 8 SUM=SUM+WEIGHT(J)*F(D) GO TO 5 8 CONTINUE SUM=SUM+WEIGHT(J)*F(ZI(J)*C+D) SUM=SUM+WEIGHT(J)*F(-ZI(J)*C+D) 5 CONTINUE C C.....MAKE INTERVAL CORRECTION AND RETURN...... GAUSS=C*SUM RETURN END DOUBLE PRECISION FUNCTION F(PSISTR) C C THE FUNCTION F COMPUTES THE THREE FUNCTIONS F1,F2 AND F3,FOR C A GIVEN SPHERICAL DISTANCE "PSISTR".THEN, C F=F1*F2*F3. C M. NAJAFI, 10/1980 IMPLICIT REAL*8(A-H,O-Z) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CZ/ ZPP(20,15),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M INDEX=INDEX+1 INFT=INFINT F=0.D0 PI=DARCOS(-1.D0) IF(PSISTR.EQ.0.D0.OR.PSISTR.EQ.PI) RETURN DR=PI/180.D0 SCALE=180.D0/PSI0 PHI=PSISTR*SCALE/DR F1=DSIN(PHI*DR) F2=P(N,M,PHI) F3=ZPP(INDEX,M+1) CALL TRAPS(0,0,1000,0,0) F=F1*F2*F3 RETURN END DOUBLE PRECISION FUNCTION P(NN,MM,TETA) C C THE FUNCTION P TAKES OUT A VALUE OF ASSOCIATED LEGENDRE C FUNCTION CORRESPONDING TO THE GIVEN ANGLE "TETA",DEGREE "NN" C AND ORDER "MM" FROM THE TABLE (ARRAY) OF THE LEGENDRE FUNCTIONS C PNM(SEE THE MAIN PROGRAM) BY A LINEAR INTERPOLATION. C M. NAJAFI, 10/1980 IMPLICIT REAL*8(A-H,O-Z) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM ANG=TETA K=(NN+1)*(NN+2)/2.D0-NN+MM IF(TETA.GT.90.D0) GO TO 99 IR=ANG/PSINT+1 IF(IR.EQ.IPNM) GO TO 88 DANG=ANG-(IR-1)*PSINT P=PNM(IR,K)+(PNM(IR+1,K)-PNM(IR,K))*DANG/PSINT GO TO 77 88 P=PNM(IR,K) 77 RETURN 99 ANG=180.D0-TETA IR=ANG/PSINT+1 DANG=ANG-(IR-1)*PSINT P=PNM(IR,K)+(PNM(IR+1,K)-PNM(IR,K))*DANG/PSINT P=P*(-1)**(NN+MM) RETURN END DOUBLE PRECISION FUNCTION H(ID,IO) C C THE FUNCTION H COMPUTES A MULTIPLIER,TO NORMALIZE AN C ASSOCIATED LEGENDRE FUNCTION OF DEGREE "ID" AND ORDER "IO". C THAT IS, _ C P =H*P C M. NAJAFI, 10/1980 IMPLICIT REAL*8(A-H,O-Z) IF(IO.NE.0) GO TO 5 H=DSQRT(2.D0*ID+1.D0) RETURN 5 NPM=ID+IO XX=2*(2.D0*ID+1.D0) ITERMS=2*IO DO 10 I=1,ITERMS XX=XX/NPM 10 NPM=NPM-1 H=DSQRT(XX) RETURN END SUBROUTINE LEGDRE(INORM,PHI,PN,NROW,NP1) C C C SUBROUTINE LEGDRE COMPUTES THE ASSOCIATED LEGENDRE FUNCTIONS C UP TO AND INCLUDING DEGREE AND ORDER N TO BE SPECIFIED. THE C DIMENSION OF PN IS PN(N+1,N+1). THE ASSOCIATED LEGENDRE POLY- C NOMIAL OF DEGREE A AND ORDER B IS STORED IN PN(A+1,B+1) IF A.NE.B C (ZONAL AND TESSERAL) AND IN PN(A+1,A+1) IF A = B (SECTORIAL) C C INORM = FLAG FOR NORMALIZATION: 1 = YES; 0 = NO C PHI = LATITUDE IN DEGREES C NROW = ROW DIMENSION OF PN IN THE CALLING PROGRAM C NP1= N + 1 = THE DIMENSION OF PN IN THE SUBROUTINE. SEE N BELOW C N = DEGREE OF LEGENDRE POLYNOMIAL C C C IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION PN(NROW,NP1) DATA DR/0.017453293D0/ N=NP1-1 N1=NP1 PHI=PHI*DR X =DSIN(PHI) DO 9 I=1,NROW DO 9 J=1,NP1 9 PN(I,J) = 0.0D0 PN(1,1) = 1.0D0 PN(2,1) = X PN(2,2) = 1.0D0 DO 15 I = 2,N I1 = I + 1 DO 10 J = 1,I1 IF(J .NE. 1) GO TO 11 PN(I1,J) = (1./I)*(X*(2.*I-1)*PN(I1-1,J) - (I-1)*PN(I1-2,J)) GO TO 10 11 IF(J.NE.2 .AND. I.NE.2) GO TO 12 PN(I1,J) = PN(I1-2,J) + (2.*I-1)*PN(I1-1,J-1) GO TO 10 12 IF(J.NE.(I1-1).AND.J.NE.I1) GO TO 13 PN(I1,J) = (2.*I-1)*PN(I1-1,J-1) GO TO 10 13 PN(I1,J) = PN(I1-2,J) + (2.*I-1)*PN(I1-1,J-1) 10 CONTINUE 15 CONTINUE DO 16 I = 1,N I1 = I + 1 DO 16 J = 1,I J1 = J + 1 COSPHI = DCOS(PHI) IF(COSPHI.EQ.0.D0) GO TO 20 GO TO 21 20 PN(I1,J1) = 0.D0 GO TO 16 21 PN(I1,J1) = (COSPHI**J)*PN(I1,J1) 16 CONTINUE IF(INORM.EQ.0) GO TO 14 N1 = N + 1 DO 17 I = 1,N1 17 PN(I,1) = DSQRT(2.D0*(I-1) + 1)*PN(I,1) DO 18 I = 1,N I1 = I + 1 DO 18 J = 1,I J1 = J + 1 NMM = I - J NPM = I + J FTERM = FACT(NMM) / FACT(NPM) NPOWER = NMM - NPM IF(MOD(NPOWER,2).EQ.0) GO TO 502 NPOWER = NPOWER + 1 FTERM = FTERM / 10D0 502 FTERM = DSQRT(FTERM) NPOWER = NPOWER / 2 FTERM = FTERM*10.**NPOWER 18 PN(I1,J1) = DSQRT(2.D0*(2.D0*I+1.D0))*FTERM*PN(I1,J1) 14 RETURN END FUNCTION FACT(N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION FA(59),FB(59),FMAN(118),NEXP(118) EQUIVALENCE (FA(1),FMAN(1)),(FB(1),FMAN(60)) DATA (FA(I),I=1,15)/ 0.1D0,0.2D0,0.6D0,0.24D0,0.12D0,0.72D0, #0.504D0,0.4032D0,0.36288D0,0.36288D0,0.399168D0,0.4790016D0, #0.62270208D0,0.871782912D0,0.1307674368D0/ DATA (FA(I),I=16,30)/ 0.20922789888 , 0.355687428096 , #0.6402373705728 , 0.12164510040883 , 0.24329020081766 , #0.51090942171709 , 0.11240007277776 , 0.25852016738885 , #0.62044840173324 , 0.15511210043331 , 0.40329146112661 , #0.10888869450418 , 0.30488834461171 , 0.88417619937397 , #0.26525285981219/ DATA (FA(I),I=31,45)/ 0.82228386541779 , 0.26313083693369 , #0.86833176188119 , 0.29523279903960 , 0.10333147966386 , #0.37199332678990 , 0.13763753091226 , 0.52302261746660 , #0.20397882081197 , 0.81591528324790 , 0.33452526613164 , #0.14050061177529 , 0.60415263063374 , 0.26582715747884 , #0.11962222086548/ DATA (FA(I),I=46,59)/ 0.55026221598121 , 0.25862324151117 , #0.12413915592536 , 0.60828186403427 , 0.30414093201713 , #0.15511187532874 , 0.80658175170944 , 0.42748832840600 , #0.23084369733924 , 0.12696403353658 , 0.71099858780486 , #0.40526919504877 , 0.23505613312829 , 0.13868311854569 / C DATA (FB(I),I=1,15)/ 0.83209871127414 , 0.50758021387722 , #0.31469973260388 , 0.19826083154044 , 0.12688693218588 , #0.82476505920825 , 0.54434493907744 , 0.36471110918189 , #0.24800355424368 , 0.17112245242814 , 0.11978571669970 , #0.85047858856786 , 0.61234458376886 , 0.44701154615127 , #0.33078854415194/ DATA (FB(I),I=16,30)/ 0.24809140811395 , 0.18854947016660 , #0.14518309202829 , 0.11324281178206 , 0.89461821307830 , #0.71569457046264 , 0.57971260207474 , 0.47536433370128 , #0.39455239697206 , 0.33142401345653 , 0.28171041143805 , #0.24227095383673 , 0.21077572983795 , 0.18548264225740 , #0.16507955160908/ DATA (FB(I),I=31,45)/ 0.14857159644818 , 0.13520015276784 , #0.12438414054641 , 0.11567725070816 , 0.10873661566567 , #0.10329978488239 , 0.99167793487095 , 0.96192759682482 , #0.94268904488832 , 0.93326215443944 , 0.93326215443944 , #0.94259477598383 , 0.96144667150351 , 0.99029007164861 , #0.10299016745146/ DATA (FB(I),I=46,59)/ 0.10813967582403 , 0.11462805637347 , #0.12265202031961 , 0.13246418194518 , 0.14438595832025 , #0.15882455415227 , 0.17629525510902 , 0.19745068572211 , #0.22311927486598 , 0.25435597334722 , 0.29250936934930 , #0.33931086844519 , 0.39699371608087 , 0.46845258497543 / C DATA (NEXP(I),I=1,30)/1,1,1,2,3,3,4,5,6,7,8,9,10,11,13,14,15,16, #17,19,20,22,23,24,26,27,29,30,31,33/ DATA (NEXP(I),I=31,73)/34,36,37,39,41,42,44,45,47,48,50,52,53,55, #57,58,60,62,63,65,67,68,70,72,74,75,77,79,81,82,84,86,88,90,91,93, #95,97,99,101,102,104,106/ DATA (NEXP(I),I=74,118) /108,110,112,114,116,117,119,121,123,125, #127,129,131,133,135,137,139,141,143,145,147,148,150,152,154,156, #158,160,162,164,167,169,171,173,175,177,179,181,183,185,187,189, #191,193,195/ C IF(N.LT.0) GO TO 51 IF(N.EQ.0) N = 1 IF(N.GT.118) GO TO 52 FACT = FMAN(N) N = NEXP(N) RETURN 51 WRITE(6,301) N 301 FORMAT(/////,'FACTORIAL FUNCTION HAS BEEN GIVEN A NEGATIVE ARGUMEN #T: N=',I6) STOP 52 WRITE(6,302) N 302 FORMAT(/////,'ARGUMENT OF THE FACTORIAL FUNCTION IS GREATER THAN #118: N=',I10) STOP END SUBROUTINE GRAPH(X,Y,NNN,XNAME,YNAME) INTEGER NNN REAL*8 X(NNN),Y(NNN),PI,RTD CHARACTER*8 XNAME,YNAME CHARACTER*1 LINE(51) REAL*8 YMIN,YMAX,YRANGE INTEGER COLUMN,XAXIS,YPRINT,I PI=DARCOS(-1.D0) RTD=180.D0/PI C FIND RANGE OF Y TO BE PLOTTED CALL MINMAX(Y,NNN,YMIN,YMAX) IF(DABS(YMIN).LE.1.D-15.AND.DABS(YMAX).LE.1.D-15) THEN DO WRITE(6,11) 'KERNEL IS ZERO' RETURN END IF 11 FORMAT(' ',//,75X,A14) YRANGE=YMAX-YMIN C CLEAR LINE TO BLANKS DO 10 COLUMN=1,51 10 LINE(COLUMN)=' ' C PLACE X-AXIS MARK IN LINE XAXIS=50*(0-YMIN)/YRANGE+1.5 IF(XAXIS.GE.1.AND.XAXIS.LE.51) THEN DO LINE(XAXIS)='I' END IF C LABEL THE GRAPH WRITE(6,15) 'GRAPH OF',YNAME,'VERSUS',XNAME 15 FORMAT('1',A8,A8,5X,A6,A8) WRITE(6,*) 'MINMUM OF',YNAME,YMIN WRITE(6,*) 'MAXIMUM OF',YNAME,YMAX,' S.D.(DEGRE +ES)',' KERNEL' C PREPARE AND PRINT LINES OF GRAPH DO 30 I=1,NNN YPRINT=50*(Y(I)-YMIN)/YRANGE+1.5 LINE(YPRINT)='*' WRITE(6,20) LINE,X(I),Y(I) 20 FORMAT('0',5X,51A1,10X,F8.2,3X,F20.12) IF(YPRINT.EQ.XAXIS) THEN DO LINE(YPRINT)='I' ELSE DO LINE(YPRINT)=' ' END IF 30 CONTINUE RETURN END SUBROUTINE MINMAX(NUMBER,NN,MIN,MAX) INTEGER NN REAL*8 NUMBER(NN),MIN,MAX MIN=NUMBER(1) MAX=MIN DO 10 I=2,NN IF(NUMBER(I).LT.MIN) MIN=NUMBER(I) IF(NUMBER(I).GT.MAX) MAX=NUMBER(I) 10 CONTINUE RETURN END