- Email: [email protected]
Computation of gravity effects of two- and three-dimensional bodies with radial symmetry
Compare & Gtoscienc¢~, Vol. 5, pp. 313-323 © Pelllam~ Press Ltd., 1979. Printed in Great Britain
009~-3~1201.-~313/$02.~
COMPUTATION OF GRAVITY EFFECTS OF TWO- AND THREE-DIMENSIONAL BODIES WITH RADIAL SYMMETRY M. M. PAm'~ and J. GOV~D~JAN Department of Physics and Ceno'efor Chemical Physics, Universityof Western Ontario, London, Ontario, Canada N6A 5B7 (Received 15 December 1978)
Almract--A computerprogram in the form of a subroutinewhichcomputesthe gravityeffectsof bodies with radial symmetry is descn'bed and presented. The program can be utilized in the situation of a circular disc, an infinite cylinder, a finite cylinder and for any arbitrary object which can be approximated into a set of circular discs about a common axis. The program has been applied to known situations and the results obtained are in excellent agreement with the previously known and reported ones. Key Words: Geophysics,Gravity surveys, Modeling. ~ntoouc-noN Intrepretation of gravity survey results has come a long way since the pioneering work of Nettleton almost four decades ago 0slettleton, 1940, 1942). With the advent of high-speed digital computers, a comparison of the gravity fields due to alternative mass distributions can he carried out efficiently. These model fields can be compared with each other as well as with the experimentally observed fields. Recently, some computer programs (Milsom and Worthington, 1977) have been published for computing the gravity effects of bodies of uniform density "and having a polygonal cross section, However, in many situations, such as salt domes, plugs of intrusive rocks, mine shafts, etc. it may be more attractive to exploit any radial symmetry that may be present in the object whose field is to be determined. If one can conceive these objects to be made up of circular discs, finite cylinder, infinite cylinder (i.e. one in which the length is much larger than the radius), hemispherical, conical, or a composite system of two or more such shapes, then the program described in this paper can be utilized effectively, The program described here can not only perform the aforementioned calculations efficiently, but also can be adapted easily to other special situations, such as when the system is not homogeneous but the density changes along the axis. The main step is the calculation of g,, the vertical component of the gravitational attraction of a circular disc at any arbitrary point. This is done from the formula described in the following section, and requires the catculation of Bessel function; the exponential function and one-dimensional integration, The programs have been written in FORTRAN IV and have been tested for a variety of situations at the University of Western Ontario, on a CYBER 73 computer, For situations, where a schematic pictorial representation of the results may be important, we have included
a subroutine GRPLT, which can draw sketches of upto four functions on the line printer. In situations where a full three-dimensional relief pattern is required, a call to the three-dimensional plotting can be made from the main program after a calculation of g, has been made. (A three-dimensionai plotting pacluge such as the one descrihed in Communications of the ACM vol. 15, p. 100, 1972may be used.) bU~Tm3AtA~C~ The vertical component of the gravitational attraction g, of a circular disc can be obtained in a number of different manners. For an arbitrary density distn'bution over the disc, one can start with the eigenfunctions of the corresponding integral operator, following the ° analogy to the corresponding electrostatic problem (Wolfe, 1971). The eigenfunctions for this situation are the Associated Legendre Polynomials. Parasnis (1961) obtained expressions for g,, for a uniform circular disc as well as an infinite cylinder in terms of an infinite series of Legendre Polynomiais. In his actual calculation, the expansion included only terms up to the fourth order. There have been a number of efforts to obtain closed form expressions (Nabighiyan, 1962; Rao and Radhakrishnamurthy, 1966; Singh, 1977a, 197To) in terms of well-known functions that have been tabulated. From the point of view of programming convenience, it is preferable to start with the Poisson's equation for the gravitational potential and numerically perform the integration in the representation of the field as an integral. As shown in Figure 1, let (r, z) be the coordinate of an arbitrary point measured from the center of the disc as origin, and the plane of the disc as the z = 0 plane. If a is the radius of the disc, p its density and t its thickness then the potential u(r, z) satisfies the equation a2u/a~ 2 "~"11 OU/OZ "~"O2U/Ot 2 = 4~Gotg(z)O(a - r) (1)
1"Permanent address: Department of Physics, Indian Institute of Technology,Kanpur, India.
where G is the gravitational constant, ~ and O are the
313
314
M.M. P~dcrand J. GOVIHDARAJAH Dirac delta function and the O function of Heaviside, respectively. Following Bateman (1959), p. 409) we have z~
z-o
i
[
;
J
and u(r, z) = -2~rGpta f ; Jl(ap)Jo(rp) e-PZlp dp,
(2)
gz = Ou/Oz = 2~Gpta f o Jl(ap)Jo(rp) e -~ dp.
(3)
The integrals in equations (2) and (3) can be represented in terms of elliptic functions, but in the program e presented in the paper, we directly integrate this numerically after generating the required Bessel func{ ~ z tions. These are calculated from the expansion formula Zo as given in Abramowitz and Stegun (1964). The upper ~ ~ z'° limitof the integral in the equations (2) and (3) does not create any practical problems for values of z, that are not extremely small. One can under these circumstances set the upper limit to 15.0/z and maintain good accuracy. In most practical geophysical prospecting one is not c ~ interested in very small or very large values of r and z, 0 but in the "small to moderate" distances from the origin Figure 1. Coordinate System for various shapes; A, disc; B, of any gravitational anomaly. Only for very small values cylinder(finiteor inf~te); C, themispere;D, cone. of z, say when z/a is less than 0.001, there is a practical .................
NORHAL I ZE0
GZ
0.0 0.0
1.0
+÷++4-÷,l.4.÷4.÷÷4-+++÷++4.~.4.+4.÷÷÷+÷+++÷ + + + 4-
÷++~+÷~÷÷~÷÷+++~÷÷÷÷÷
++++÷÷~÷÷~÷+~+÷~÷÷÷~+~÷÷÷~++~ A + + 4. BA +
÷ 4.
÷ +
÷
+ + + 4. + + 4. 4,. ,~
OC
D C
+ 4, +
0C
8CA
RIA
4. 4. ÷ B
S
C
4.
+ +
C
D
O
+ + + + + + ÷ 4.
A 4. 4. +
~8
A + + +
+
4. A 4. +
4. 4. A B 4. 4"
A ~+ 4. A 4. + 4. +
+ +
A 4. +
A 4. ÷
D
S
8
B
B
S
C
C
D
0
+ ÷ 4" + + ÷ + + +
C
D
÷
+ + + C O
+
Figure2. Line-printerplotof normalizedgz forcirculardiscas functionof r/a (A, B, C, D representincreasingvalues of z/a).
Computation of gravity effects of two- and three-dimensional bodies problem in evaluating the integral because of the Bessel functions exhibiting oscillatory behavior for upto very large values of p. There are a number of well-known methods for such oscillating integrands, such as the Fiion quadrature, the Euler transformation or the transformarion to a new set of variables, However, because one can change the number of integration mesh points, one need not have to take recourse to any of the special methods mentioned previously. What the program does is to perform a more accurate integration, when z is smaller than a certain specified value. This arrangement is required only for the situation of a disc, because for the cylinder, we evaluate the expression (2) and further due the presence of p in the denominator the oscillations are damped anyway. For the situation of an infinite cylinder, we choose the top of the cylinder to be at z = 0, and let the cylinder extend to z = -®. For a cylinder of finite length 1, the integral in equation (2) has to be evaluated twice, once with z equal to the z-coordinate of the point of interest, followed by a z ~ 1 + z, the coordinate of the point of interest. The difference of the two yields gz at (r, z) for a finite cylinder. For any other three-dimensional object with radial symmetry the top surface may be chosen as the z = 0 plane and its length divided into a number of segments, to each of which equation (3) is applicable and the net effect determined by an integration along z. US~GTHEPaOGmLM To calculate the integrals in equations (2) and (3) all one has to do is to write a main program as shown in the test example and make a call to the subroutine AXIAL, with the following parameters: CALL AXIAL(R, Z, RZ,, RL, RS) where the parameters in the call have the following meaning. R and Z are the r, z coordinates of the point of interest, RZ is the radius of the disc, or of the cross section of the cylinder, RL is the length of the cylinder for the finite cylinder case and RS is the output from this subroutine and represents the following dimensionless form of the integral in equation (3), fo RS =
Yl(t)Jdxt)e-" dt
(4)
where x = r/a and y = zla. In all practical applications, one would be interested in calculating gz at a large number of points for which experimental observations are planned or already available. Hence a short main program somewhat along the lines of the test program which we present in the Appendix would have to be written to define the points of interest. This main program must have the common declaration (~OMMON X, I/, ISHAPE and ISHAPE must be defined before reference is made
315
to the subroutine AXIAL. This parameter ISHAPE must be assigned a value of 1, 2, or 3 depending on whether calculations are to be made for a circular disc, an infinite cylinder, or a finite cylinder respectively. The output RS returned by the subroutine AXIAL may be multiplied by the appropriate units in the main program to give g~ at the required point. The subroutine AXIAL calls two functions PINT and FINK. The function FINK performs the task of calcularing the integrand in equation (4) for arbitrary values of T for a given X and Y. The X and Y are referenced through the blank common mentioned previously. The function RINT to which a call is made by the statement RS = RINT (FINK, A, B, EPS) integrates the function from A to B by an adaptive Simpson's quadrature. EPS is an error parameter, and can be set to 0.0001. The lower limit for the integral A is always set to 0, and B to 15.0/¥, for most situations. For very small values of Y, the integration is bypassed and the result is set to the analytically known result for Y = 0. tmULTS~OM nOGRAM The main aim of this article is to provide a program which calculates gz at the point (r, z) for objects with radial symmetry, thus reducing the labor of looking into tabulated functions and the interpolations. We therefore have used the program to obtain illustrative results for special situations, for which exact results have been reported earlier in the literature. As a first application, we have calculated the gz at various (r, z) for a circular disc. Our test program results reproduce exactly the results as presented by Singh (1977a, table 1). We did calculations at a number of other points for display from GRPLT with a line printer as well as with the CALCOMP plotter. Our results are somewhat at variance with those reported by Parasnis (Parasnis, 1961), at times by about 5 percent. If anything, this implies that his truncation at the fourth term yields reasonable results. Our direct integration technique is more versatile, convenient and accurate than an expansion in terms of Legendre functions. We also have used the program for the situation of an infinite cylinder, and finite cylinder as well as for an hemisphere and a cone. For points on the axis, the results agree closely with those from the analytic ones for these situations. For a compilation of the appropriate expressions the reader is referred to the review article by Talwani (1973). ADArrA~oNsANDEXTENSIONS The program presented here is applicable as such for a large number of situations of practical intero, t. With a suitably written program segment it may be used for any composite shape which can be built up of cylinders, cones, and spheres. By allowing for different densities of the discs into which the three-dimensional body may be divided, the program may be extended to take account of the axial variations in the density.
316
M.M. PANTand J. G O v m D ~
~,ud£NCES AMamowitz, M., and Stegun, I. A., 1964, Handbook of mathcmatical functions: Dover Publ., Inc., New York, p. 369-370. Bateman, H., 1959, Partial differential equations of mathematical physics: Cambridge Univ. Press, Cambridge, 469p. Milsom, J., and Worthington, G. A., 1977, Computer pmgralns for rapid computation of gravity effects of two-dimensional and three-dimensional bodies: Computers & Geosciences, v. 3, no. 2, p. 269-281. Nabiighian, M. N., 1962, The gravitational attraction of a right vertical circular cyfinder at points external to it: Geof. pura e Appl. v. 53, p. 45-51. Nettleton, L. L., 1940, Geophysical prospecting for oil: McGraw Hill Book Co., New York, 444 p. Nettleton, L. L., 1942, Gravity and magnetic calculations: Oeophysics, v. 7, no. 3, p. 293--310. Parasnis, D. S., 1961, Exact expressions for the gravitational
attraction of a circular lamina at all points of space and of a fight circular vertical cylinder at points external to it: Oeophys. Prosp., v. 9, p. 382-398. Rao, B. S. R., and Radhakrishnamurty, 1966, An expression for the gravitational attraction of circular plate: Indian Jour. Pure Appi. PhD., v. 4, p. 276--278. Singh, S. K., 1977, Gravitational attraction of a circular disc: Geophysics, v..42, no. 1, p. 111-113. Singh, S. K., 1977a, Gravitational attraction of a vertical right circular cylinder! Geophys. Jour. Roy. Astr. Sc~., v. 50, p. 243-246. Talwani, M. 1973, Computer usages in the computation of gravity anomalies,/n Methods of computational physics, v. 13: ed. by B.A. Bolt, Academic Press, New York, p. 343--389. Wolfe, P., 1971, Eigentunctions of the integral equation for the potential of the charged disk: Jour. Math. Phys., v. 12, p. 1215-1222.
AEPENI)C(
PROGRAM TEST(INPUT,OUTPUT,TAPES=INPUT,TAPE6=OUTPUT,PUNCH) COMMON X,Y, ISHAPE DIMENSION S X C l l 0 ) , S Y C l l 0 ) , S F B C I 1 0 , 1 1 0 ) , S P H C I 0 0 ) , D V ( 5 0 , ~ ) C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
THIS IS A TEST PROGRAM WHICH IS DESIGNED TO SHOW HOW THE SUBROUTINE AXIAL AND THE RELATED ROUTINES ARE TO BE ACCESSED IN ORDER TO CALCULATE THE Z COMPONENT OF THE GRAVITATIONAL F I E L D . THERE ARE FOUR SECTIONS OF THIS PROGRAM WHICH DEALS WITH THE CALCULATIONS FOR OBdECTS SUCH AS A D I S C , I N F I N I T E CYLINDER,HEMISPHERE AND A RIGHT CIRCULAR CONE RESPECTIVELY. THE VARIABLE ISHAPE IN THE BLANK COMMON IS TO BE SET BY THE USER IN THE MAIN PROGRAM ACCORDING TO THE FOLLOWING DESCRIPTION. ISMAPE
SX SY NX NY
= 1 2 3
FOR CIRCULAR DISC I N F I N I T E CYLINDER F I N I T E CYLINDER
ARRAY OF R COORDINATES ARRAY OF Z COORDINATES NUMBER OF ELEMENTS IN THE SX ARRAY NUMBER OF ELEMENTS IN THE SY ARRAY
NDISC NUMBER OF DISCS INTO WHICH THE OBdECT IS TO BE D I V I D E D . NOTE THAT THIS IS NOT NECESSARY IN THE CASE OF A DISC OR A CYLINDER. IN THE TEST EXAMPLE DESCRIBED HERE NDISC HAS BEEN ARBITRARILY SET TO l O 0 . DSX AND DSY DE4OTE THE INTERVAL BETWEEN TWO SUCCESSIVE VALUES OF SX AND SY RESPECTIVELY. IPLT .GT. .LE.
CONTROLS THE PLOT ON THE L I N E 0 0
PRINTER
PLOT IS PRODUCED PLOT 1S OMITTED
FOR THE CASE OF A RIGHT CIRCULAR CONE THE FOLLOWING THREE PARAMETERS MUST BE GIVEN AS INPUT. A BETA ZER
RADIUS OF THE BASE OF THE CONE HALF THE ANGLE OF THE CONE I N DEGREES DENOTES Z0 IN THE FIGURE AND CORRESPONDS TO THE HEIGHT OF THE APEX OF THE CONE ABOVE THE Z=0 PLANE
ITEST CONTROL VARIABLE FOR THIS TO PERFORM S P E C I F I C TASKS =l 2
C
3
C C C C C
~
TEST PROGRAM ONLY
DISC I N F I N I T E CYLINDER HEMISPHERE CONE
FOR THE CASE OF A HEMISPHERE, THE RADIUS OF THE HEMISPHERE IS TO BE GIVEN AS THE INPUT. ( REPRESENTED BY THE VARIABLE RHS )
Compu~tionofgmvitye~ectsoftwo-~d ~e~imension~bodies
317
C C C C C C C C
THE OUTPUT FROM THE SUBROUTINE AXIAL IS CONTAINED IN THE VARAIBLE RS. THE VARIABLE SFB IS A MATRIX WHICH IS USED TO COPY THE RESULT AS A FUNCTION OF RZERO AND Z RESPECTIVELY. WRITE(6,10) 10 FORMAT(IHI) 20 CONTINUE READ(5,30) ITEST,NX,NY,IPLT,ND1SCpDSX,DSY 30 F O R M A T ( S I S , 2 F I 0 . 5 )
C C C C C
THE FOLLOWING STATEMENT CHECKS FOR THE END OF F I L E ON INPUT AND IS APPLICABLE TO ONLY CONTROL DATA MACHINES. FOR OTHER COMPUTERS CHECK THE FORTRAN LANGUAGE MANUAL.
IF(EOF(5)) 3 1 0 , ~ 0 W0 CONTINUE WRITE(6,SD) ITEST,NX,NY,IPLT,NDISC,DSX,DSY 50 FORMAT(IHO,SHITEST,IX,I),2X,3HNX ,I5,1X,)HNY , I 5 , 1 X , S H 1 P L T x 6HNDISC ,15,1X,3flDSX,FI0.5,1X,3HDSY,FI0.5) C C C
GENERATE
60 70
80
90 100
C C C C C
,15,1X,
SX AND SY ARRAYS
IF(NDISC.LE.0) NDISC =100 SX(1)=0.0 SY(1)=0.0 DO 60 I=2,NX SX(1)=SX(I-I)+OSX CONTINUE DO 70 I =2,NY SY(1)=SY(I-I)+DSY CONTINUE GO TO (80, I~0, 160, 210) ITEST CONTINUE ISHAPE = 1 DO 100 I = 1,NX DO 90 d = 1,NY R = SX(1) Z = SY(J) CALL A X I A L ( R , Z t l . 0 , 1 . 0 , R S ) SFB(I,J) = RS CONTINUE CONTINUE IF(.NOT.(IPLT.GT.0)) GO TO 280 WRITE(6,10) M = 50 OUTPUT PLOT IS WRITTEN ON UNIT 6 . IF THE PLOT IS DESIRED ON ANY OTHER UNIT , CHANGE NC TO THE APPROPRIATE VALUE
NC:6 N=k DO 130 d = l , N DO 120 I = I , M D V ( l , d ) = SFB(I,J)ISFB(I,d) 120 CONTINUE 130 CONTINUE CALL G R P L T ( N C , S X , D V ~ M p N ) GO TO 280 1W0 CONTINUE C C C
INFINITE
CYLINDER
DO 155 I = I , N X OO 150 J = I , N Y R = SY(1) Z = SY(d) ISHAPE = 2 CALL A X I A L ( R , Z , I . 0 , 1 . 0 , R S ) SFB(I,J)=RS 150 CONTINUE 1SS CONTINUE GO TO 280 160 CONTINUE C
¢
318
M.M. PANTand J. GOVINDAP~JAN C C C
HEMISPHERE READ THE RADIUS OF THE HEMISPHERE 165
170
180 190 200 210 C C C
READ(5,165) RHS FORMAT(F10.5) NDM = NDISC - 1 DO 200 KK = I,NX R = SX(KK) DO 190 K = 1,NY XF = SY(K) DX=RHS/FLOAT(NDISC) DO 170 I = I , N D I S C XX=(l-l)xOX RZ=SQRTCXXX(2.OXRHS-XX)) Z=XF+RHS-XX ISHAPE = 1 CALL A X I A L ( R , Z , R Z , I . 0 , R S ) SPH(1) = RS CONTINUE ZPP = 0 . 5 ; : ( S P H ( 1 ) + SPH(NDISC))XDX OO 180 I=IpNOM ZPP=ZPP+SPHCI~XDX CONTINUE SFB(KK,K) = ZPP CONTINUE CONTINUE GO TO 280 CONTINUE RIGHT CIRCULAR CONE
READ(S,220) A, BETA,ZER 2"20 FORMAT(3F10.5) WR1TEC6,230) A,BETA,ZER 230 FORMAT(IH0,1OX,3H A = , F 6 . 2 , 2 X , 6 H C C C
BETA=,2X,FI0.2,2X,SH ZER=,FL0.3,/)
CONVERTBETA FROM DEGREES TO RADIANS PI=q.0XATAN(I,0) BETA=BETA;PI/180.0 NDM = NDISC - 1 CT=I.0/TANCBETA) DX=(AXCT-ZER)/NOISC DO 270 1 = I,NX DO 260 J = 1,NY R : SX(l) DO'2~0 K = 1,NDISC RZ =( ZER + K~DX)XTAN(BETA) Z = SY(J) + KXOX ISHAPE = 1 CALL A X I A L ( R , Z , R Z , I . 0 , R S ) 2~0 SPH(K)=RS SUM=0.5~DXX(SPH(1)+SPH(NDISC)) DO 250 M=2,NDM SUM=SUM+OXXSpH(M)
250 CONTINUE S F B ( 1 , J ) = SUM 260 CONTINUE 270 CONTINUE 280 CONTINUE WRITE(6,10) DO 285 I=I,-NX WRITE(6,290) C S F B C I , J ) , J = I , N Y ) 285 CONTINUE 290 FORMATCIHO,CSF15.6,/)) PUNCH 3 0 0 , ( ( S F B ( 1 , J ) , I = I , N X ) , J = I , N Y ) 300 FORMAT(8F10.6) GO TO 20 310 CONTINUE STOP END SUBROUTINE A X I A L ( R , Z , R Z , R L , R E S ) COMMON X , Y , ISHAPE EXTERNAL FINK C C C C
THIS SUBROUTINE IS THE KEY SUBROUTINE WHICH ACTS AS THE INTERFACE BETWEEN THE USER DEFINED MAIN PROGRAM AND THOSE ROUTINES WHICH PERFORM THE ACTUAL CALCULATION
Computation of gravity effects of two- and three-dimensionalbodies C C C C C C C C C C C C C C
319
INPUT . . . . . . . . . . . R AND Z ARE RESPECTIVELY THE COORDINATES OF THE POINT AT WHICH THE F I E L D IS TO BE CALCULATED RZ
THE RADIUS OF THE DISC OR THE CYLINDER
RL THE LENGTH OF THE F I N I T E
CYLINDER
OUTPUT . . . . . . . . . RES
VARIABLE
WHICH CONTAINS THE RESULT
YZ = 0 . 0 0 1 X=R/RZ Y=Z/RZ A=0.0 B = 20.0/Y EPS=0.0001 IF(Y.LT.YZ) GO TO 10 IF (Y . L T . 0 . 0 1 5 ) EPS = 0 . 0 0 0 0 1 RES=Ri N T ( F I N K , A, B, EPS) IF(ISHAPE.LE.2) GO TO 20
Y=(Z+RL)/RZ A = 0.0 B = 20.0/Y EPS = 0 . 0 0 0 1 RESA=RI NT(F I N K , A , B , E P S ) RES=RES - RESA GO TO 20 10 CONTINUE IF(X.LT.1.0) RES = 1 . 0 IF(X.EQ. 1.0) RES=0.5 IF(X.GT.1.0) RES=0.0 20 CONTINUE RETURN END
FUNCTION R I N T ( F , A , B t E P S ) DIMENSION D X C 3 0 ) , E P S P C 3 0 ) , X 2 C 3 0 ) , X 3 C 3 0 ) , F 2 C 3 0 ) , F 3 ( 3 0 ) , F ~ ( 3 0 ) ,
=
FMP(~0),FBP(~0),EST2(~O),EST3(~O),PVAL(~0,3) INTEGER RTRNC30)
C
C C C
THIS FUNCTION SUBPROGRAM PERFORMS THE INTEGRATION BY A NONRECURSIVE ADAPTIVE SINP$ONS INTEGRATION FORNULA. THE FORTRAN CODE IS BASED ON THE ALGORITHM 182
C C
PUBLISHED
IN THE COMMUNICATIONS OF ACM, PAGE 2~Np
196k.
LVL=0 ABSA=I.0 EST=I.0 DA=B-A FA=F(A) FM=~.0xFCCA+B)/2.0) FB=F(B) 10 CONTINUE LVL=LVL+I DX(LVL)=OA/3.0 SX=DX(LVL)/6.0 FI=~.0XF(A+DX(LVL)/2.0) X2(LVL)=A+OX(LVL) F2(LVL)=F(X2(LVL)) X3CLVL)=X2CLVL)+DXCLV L) F3(LVL)=F(X3(LVL)) EPSP(LVL)=EPS Fk(LVL)=h.0;:F(X3(LVL)+OX(LVL)/2.0) FMP(LVL)=FM ESTI=(FA+FI+F2(LVL));SX FBP~LVL)=FB EST2(LVL):(F2(LVL)+F3(LVL)+FM)XSX EST3(LVL)=(F3(LVL)+F~(LVL)+FB)XSX SUM=ESTI+EST2(LVL)+EST3(LVL) ABSA=ABSA-ABS(EST)+ABS(EST1)+ABS(EST2(LVL))+ABS(EST3(LVL))
IF(.NOT.((ABS(EST-SUM).LE.EPSP(LVL)XABSA).AND.(EST.NE.1.0)) X(LVL.GE.30))GO 20 CONTINUE LVL=LVL-I
TO 30
.OR.
3~
M.M. PAN'r~dJ.OOV[NIDA~JAN
30
~0
50
60
PVAL(LVL,RTRN(LVL))=SUM GO TO ( ~ 0 , 5 0 , 6 0 ) RTRN(LVL) CONTINUE RTRNCLVL)=I DA=DX(LVL) FM=F1 FB=F2(LVL) EPS=EPSPCLVL)/1.7 EST=ESTI GO TO 10 CONTINUE RTRNCLVL)=2 DA=DX(LVL) FA=F2(LVL) FM=FHP(LVL) FB=F3(LVL) EPS=EPSP(LVL)/I.7 EST=EST2(LVL) A=X2CLVL) GO TO 10 CONTINUE RTRNCLVL)=3 DA=DX(LVL) FA=F3(LVL) FM=F~(LVL) FB=FBP(LVL} EPS=EPSP(LVL) EST=EST3CLVL) A=X3(LVL) GO TO 10 CONTINUE SUM=PVAL(LVLtl)+PVAL(LVL,2)+PVAL(LVL,3) | F C L V L . G T . I ) GO TO 20 RINT=SUM RETURN END FUNCTION FINKCT)
C C C C C C
THIS FUNCTION SUBPROGRAM CALCULATES THE BESSEL FUNCTIONS BASED ON THE FORMULAE AS GIVEN IN THE HANDBOOK OF MATHEMATICAL FUNCTIONS BY ABRAMOWITZ AND STEGUN. DIMENSION A ( 6 ) , B ( 7 ) , C ( 7 ) , D ( 6 ) , E ( 7 ) , F ( 7 ) COMMON X , Y , ISHAPE DATA A I - 2 . 2 ~ 9 9 9 9 7 , 1 . 2 6 5 6 2 0 8 , - 0 . 3 1 6 3 8 6 6 , 0 . 0 ~ k ~ 7 9 , - 0 . 0 0 3 9 ~ , 1 0.0002100/
C DATA B / 0 . 7 9 7 8 8 ~ 5 6 , - 0 . 0 0 0 0 0 0 7 7 , - 0 . 0 0 5 5 2 7 h 0 , - 0 . 0 0 0 0 9 5 1 2 ,
1
0.0013727,-0.00072805,0.0001~k761
C 1
DATA C / - 0 . 7 8 5 3 9 8 1 6 t - 0 . 0 ~ 1 6 6 3 9 7 , - 0 . 0 0 0 0 3 9 5 ~ , 0 . 0 0 2 6 2 5 7 3 , -0.0005~125,-0.000293~3,0.00013558/
1
DATA D / - 0 . 5 6 2 k 9 9 8 5 , 0 . 2 1 0 9 3 5 7 ~ , - 0 . 0 3 9 5 ~ 2 8 9 , 0 . 0 0 ; ~ 3 3 1 9 , - 0 . 0 0 0 3 1 7 6 1 , 0.00001109/
1
DATA E 10.79788~56,0.00000156,0.01659667,0.00017105,-0.002~9511, 0.00113656,-0.00020033/
1
DATA F / - 2 . 3 5 6 1 9 k k 9 , 0 . 1 2 k 9 9 6 1 2 , 0 . 0 0 0 0 5 6 5 0 , - 0 . 0 0 6 3 7 8 7 9 , 0 . 0 0 0 7 ~ 3 h B , 0.0007982~,°0.00029166/
C C C C DATA TL I I . O E - 0 7 / C C C
MAIN LOOP GENERATING BESSEL FUNCTIONS I F C T . L E . T L ) T = TL XT=X~T yT=Y~:T IF(T-3.0) 10,10,30 10 CONTINUE X3=T/3.0 X32=X3~X3 FOI=0.5 DO 20 g = l , 6 Fdl=FJ1+D(K)=X32~;K 20 CONTINUE
Computation of gravity effects of two- and three-dimensional bodies
30
~0 50 60
70 80
90 100
Fdl=FdlXT GO TO 50 CONTINUE X3=3.0/T FFI=E(1) FF2=F(I) DO ~0 K = 2 , 7 X3K=X3X=C(K-1) FFI=FFI+E(K)XX3K FF2=FF2+FCK)XX3K CONTINUE FdI=FFI=tCOS(FF2+T)/SQRT(T) CONTINUE IF(XT-3.0) 60,60,80 CONTINUE X3=XT=XT/9.0 FdZ=I.0 DO 70 K = l , 6 FJZ=FdZ+A(K)=X3=XK CONTINUE GO TO 100 CONTINUE X3=3.0/XT FFI=B(I) FF2=C(1)+XT DO 90 K = 2 , 7 X3K=X3X~(K-1) FFI=FFI+B(K)=X3K FF2=FF2+C(K)~X3K CONTINUE FdZ=FFI=COSCFF2)/SQRT(XT) CONTINUE FINK=FdZ=Fdl~EXP(-YT) IF ( I S H A P E . E Q . 2 ) F I N K = F I N K / T IF ( I S H A P E . E Q . 3 ) F I N K = F I N K / T RETURN END SUBROUTINE GRPLT(NC,X,YsM,N)
C C C C
THIS SUBROUTINE PERFORMS THE TASK OF CREATING A PLOT ON THE LINE PRINTER.
C
NC IS THE NUHBER OF THE DATA CHANNEL
C C C C C
THI$ SUBROUTINE IS BASED ON THE ALGORITHM PUBLISHED IN THE COMMUNICATIONS OF ACM, P ~ 9 2 , 1 9 7 1 AS ALGORITHM ~ 1 2 .
C C C
USUALLY UNIT 6
DIMENSION X(M),Y(H,N),LP(50),IW(IOO),PS(IO0) INTEGER S,LsPS DATA X E R R / 0 . 0 0 0 I / , Y E R R / 0 . 0 0 0 1 / DATA L / 5 0 / , 5 / 8 5 / DATA I B / 1 H / DATA MC/1H+/ DATA I C 1 / I H A / DATA I C 2 / I H B / DATA I C 3 / 1 H C / DATA IC~/1HO/ FIND XHAX AND XMIN XMAX=X(I) XMIN=X(1) DO 10 l=2,N IF(X(1).GT.XMAX) IF(X(I).LT.XMIN) 10 CONTINUE
C C C
-
XMAX=XCI) XMIN=X(I)
FIND YHAX AND YHIN
YMIN=YCI,I) YMAX=Y(1,1) DO 50 DO 20
I=1,~ d=l,N
IF(Y(I,d).GT.YMAX) YMAX=YCI,d) IF(Y(I,d).LT.YMIN) YMIN=Y(I,J) 20 CONTINUE
321
M.M. PANT and J. GOVlND~AN
322
30 CONTINUE 00 ~0 I=I,M OO 35 d = l , N Y(I,d)=T(I,d)/YMAX 35 CONTINUE 40 CONTINUE XMIN=0.0 YMIN=0.0 YMAX=I.O IFCCABS(XMAX-XMIN).LT.XERR).OR.CABS(YMAX-YMIN).LT.YERR))RETURN
PX=FLOAT(L-1)/(XMAX-XMIN)
50 60
70 80 85
PY=FLOAT($-I)I(YMAX-YMIN) DO 50 I = I , L LP(I)=0 CONTINUE KL=M+I CONTINUE KL=KL-1 IF(KL.LT.I) GO TO 70 KR=I+INT(CX(KL)-XMIN)XPX+0.5) LP(KR)=KL GO TO 60 CONTINUE WRITE(NC,80) FORMAT(1HO) WRITECNC,$5) FORMATC1H , ~ O X , I S H NORMALIZED GZ ) DO 90 d = l , S
IW(J)=MC 90 CONTINUE WRITE(NCt95) 95 F O R M A T ( l O X , 3 H O . 0 , 8 2 X , ) H I . 0 ) WRITECNC,100) ( I W ( d ) , d = I , S ) 100 FORMAT(1H , 1 0 X , 1 0 0 A 1 ) DO 250 I = I , L PS(1)=l PS(S)=I MS=S-1 DO 110 J=2,M$ PS(d)=2 110 CONTINUE IF(.NOT.(LP(I).GT.0)) GO TO 150 120 CONTINUE IF(.NOT.(N.EQ.1)) GO TO 130
INTR=I+INT(0.5+PY::(Y(LP(1),I)-YMIN)) PS(INTR)=3 GO TO 150 130 CONTINUE d=N+l 140 CONTINUE
J=d-I IF(J.LT.I) GO TO 150 INTR=I*INT(0.5+PY;:(Y(LP(1),J)-YMIN)) PSCINTR)=J+2 GO TO 14O 150 CONTINUE DO 220 J = l , S
IW(U)=IB IF(.NOT.(PS(d).EQ.1))
GO TO 170
IW(J)=MC GO TO 220 170 CONTINUE IF(.NOT.(PS(J).EQ.2))
GO TO 180
lW(d)=I8 GO TO 220 180 CONTINUE IFC.NOT.(PSCd).EQ.3))
GO TO 190
IW(d)=ICl GO TO 220 190 CONTINUE IF(.NOT.(PS(J).EQ.~))
GO TO 200
GO TO 220 200 CONTINUE IF(.NOT.(PS(d).EQ.5))
GO TO 210
IW(J)=IC2
IW(d)=IC3 GO TO 220 210 CONTINUE
IW(d)=IC~
Computation of gravity e~ects of two- and three-dimensional bodies
220 CONTINUE IFC.NOT.CI.EQ.1)) GO TO 230 WRITECNC,225) (IWCd),d=l,S) 225 FORMATC6Xt3H0.0t2X,100A1) GO TO 250 230 CONTINUE IFC.NOT.CI.EQ.25)) GO TO 2k0 WRITECNC,235) (IWCd),J=I,S) 255 FORMAT(6X,3HR/A,2X,10OA1) GO TO 250 2k0 CONTINUE WRITECNC,2qS) ( I W ( d ) e d = l , $ ) 2k5 FORIqATC1H ,10Xpl00Al) 250 CONTINUE 260 FORMAT(1H ,11) DO 270 d=l,S IW(d)=MC 270 CONTINUE WRITECNCe275) XMAX,CIWCJ),J=lt$) 275 FORRAT(6X,Fk.2,1X,10OAI) WRITE(NC,260) WRITECNC;2$0) XMINtXMAX,YMIN,YMAX 280 FORHAT(1HOpIOX,kHXMIN,FIO.k,11X,~HXMAX,F10.~,llX,kHYMIN,FlO.qp 111X,~HYMAJ(,FI0.N) RETURN END
323
009~-3~1201.-~313/$02.~
COMPUTATION OF GRAVITY EFFECTS OF TWO- AND THREE-DIMENSIONAL BODIES WITH RADIAL SYMMETRY M. M. PAm'~ and J. GOV~D~JAN Department of Physics and Ceno'efor Chemical Physics, Universityof Western Ontario, London, Ontario, Canada N6A 5B7 (Received 15 December 1978)
Almract--A computerprogram in the form of a subroutinewhichcomputesthe gravityeffectsof bodies with radial symmetry is descn'bed and presented. The program can be utilized in the situation of a circular disc, an infinite cylinder, a finite cylinder and for any arbitrary object which can be approximated into a set of circular discs about a common axis. The program has been applied to known situations and the results obtained are in excellent agreement with the previously known and reported ones. Key Words: Geophysics,Gravity surveys, Modeling. ~ntoouc-noN Intrepretation of gravity survey results has come a long way since the pioneering work of Nettleton almost four decades ago 0slettleton, 1940, 1942). With the advent of high-speed digital computers, a comparison of the gravity fields due to alternative mass distributions can he carried out efficiently. These model fields can be compared with each other as well as with the experimentally observed fields. Recently, some computer programs (Milsom and Worthington, 1977) have been published for computing the gravity effects of bodies of uniform density "and having a polygonal cross section, However, in many situations, such as salt domes, plugs of intrusive rocks, mine shafts, etc. it may be more attractive to exploit any radial symmetry that may be present in the object whose field is to be determined. If one can conceive these objects to be made up of circular discs, finite cylinder, infinite cylinder (i.e. one in which the length is much larger than the radius), hemispherical, conical, or a composite system of two or more such shapes, then the program described in this paper can be utilized effectively, The program described here can not only perform the aforementioned calculations efficiently, but also can be adapted easily to other special situations, such as when the system is not homogeneous but the density changes along the axis. The main step is the calculation of g,, the vertical component of the gravitational attraction of a circular disc at any arbitrary point. This is done from the formula described in the following section, and requires the catculation of Bessel function; the exponential function and one-dimensional integration, The programs have been written in FORTRAN IV and have been tested for a variety of situations at the University of Western Ontario, on a CYBER 73 computer, For situations, where a schematic pictorial representation of the results may be important, we have included
a subroutine GRPLT, which can draw sketches of upto four functions on the line printer. In situations where a full three-dimensional relief pattern is required, a call to the three-dimensional plotting can be made from the main program after a calculation of g, has been made. (A three-dimensionai plotting pacluge such as the one descrihed in Communications of the ACM vol. 15, p. 100, 1972may be used.) bU~Tm3AtA~C~ The vertical component of the gravitational attraction g, of a circular disc can be obtained in a number of different manners. For an arbitrary density distn'bution over the disc, one can start with the eigenfunctions of the corresponding integral operator, following the ° analogy to the corresponding electrostatic problem (Wolfe, 1971). The eigenfunctions for this situation are the Associated Legendre Polynomials. Parasnis (1961) obtained expressions for g,, for a uniform circular disc as well as an infinite cylinder in terms of an infinite series of Legendre Polynomiais. In his actual calculation, the expansion included only terms up to the fourth order. There have been a number of efforts to obtain closed form expressions (Nabighiyan, 1962; Rao and Radhakrishnamurthy, 1966; Singh, 1977a, 197To) in terms of well-known functions that have been tabulated. From the point of view of programming convenience, it is preferable to start with the Poisson's equation for the gravitational potential and numerically perform the integration in the representation of the field as an integral. As shown in Figure 1, let (r, z) be the coordinate of an arbitrary point measured from the center of the disc as origin, and the plane of the disc as the z = 0 plane. If a is the radius of the disc, p its density and t its thickness then the potential u(r, z) satisfies the equation a2u/a~ 2 "~"11 OU/OZ "~"O2U/Ot 2 = 4~Gotg(z)O(a - r) (1)
1"Permanent address: Department of Physics, Indian Institute of Technology,Kanpur, India.
where G is the gravitational constant, ~ and O are the
313
314
M.M. P~dcrand J. GOVIHDARAJAH Dirac delta function and the O function of Heaviside, respectively. Following Bateman (1959), p. 409) we have z~
z-o
i
[
;
J
and u(r, z) = -2~rGpta f ; Jl(ap)Jo(rp) e-PZlp dp,
(2)
gz = Ou/Oz = 2~Gpta f o Jl(ap)Jo(rp) e -~ dp.
(3)
The integrals in equations (2) and (3) can be represented in terms of elliptic functions, but in the program e presented in the paper, we directly integrate this numerically after generating the required Bessel func{ ~ z tions. These are calculated from the expansion formula Zo as given in Abramowitz and Stegun (1964). The upper ~ ~ z'° limitof the integral in the equations (2) and (3) does not create any practical problems for values of z, that are not extremely small. One can under these circumstances set the upper limit to 15.0/z and maintain good accuracy. In most practical geophysical prospecting one is not c ~ interested in very small or very large values of r and z, 0 but in the "small to moderate" distances from the origin Figure 1. Coordinate System for various shapes; A, disc; B, of any gravitational anomaly. Only for very small values cylinder(finiteor inf~te); C, themispere;D, cone. of z, say when z/a is less than 0.001, there is a practical .................
NORHAL I ZE0
GZ
0.0 0.0
1.0
+÷++4-÷,l.4.÷4.÷÷4-+++÷++4.~.4.+4.÷÷÷+÷+++÷ + + + 4-
÷++~+÷~÷÷~÷÷+++~÷÷÷÷÷
++++÷÷~÷÷~÷+~+÷~÷÷÷~+~÷÷÷~++~ A + + 4. BA +
÷ 4.
÷ +
÷
+ + + 4. + + 4. 4,. ,~
OC
D C
+ 4, +
0C
8CA
RIA
4. 4. ÷ B
S
C
4.
+ +
C
D
O
+ + + + + + ÷ 4.
A 4. 4. +
~8
A + + +
+
4. A 4. +
4. 4. A B 4. 4"
A ~+ 4. A 4. + 4. +
+ +
A 4. +
A 4. ÷
D
S
8
B
B
S
C
C
D
0
+ ÷ 4" + + ÷ + + +
C
D
÷
+ + + C O
+
Figure2. Line-printerplotof normalizedgz forcirculardiscas functionof r/a (A, B, C, D representincreasingvalues of z/a).
Computation of gravity effects of two- and three-dimensional bodies problem in evaluating the integral because of the Bessel functions exhibiting oscillatory behavior for upto very large values of p. There are a number of well-known methods for such oscillating integrands, such as the Fiion quadrature, the Euler transformation or the transformarion to a new set of variables, However, because one can change the number of integration mesh points, one need not have to take recourse to any of the special methods mentioned previously. What the program does is to perform a more accurate integration, when z is smaller than a certain specified value. This arrangement is required only for the situation of a disc, because for the cylinder, we evaluate the expression (2) and further due the presence of p in the denominator the oscillations are damped anyway. For the situation of an infinite cylinder, we choose the top of the cylinder to be at z = 0, and let the cylinder extend to z = -®. For a cylinder of finite length 1, the integral in equation (2) has to be evaluated twice, once with z equal to the z-coordinate of the point of interest, followed by a z ~ 1 + z, the coordinate of the point of interest. The difference of the two yields gz at (r, z) for a finite cylinder. For any other three-dimensional object with radial symmetry the top surface may be chosen as the z = 0 plane and its length divided into a number of segments, to each of which equation (3) is applicable and the net effect determined by an integration along z. US~GTHEPaOGmLM To calculate the integrals in equations (2) and (3) all one has to do is to write a main program as shown in the test example and make a call to the subroutine AXIAL, with the following parameters: CALL AXIAL(R, Z, RZ,, RL, RS) where the parameters in the call have the following meaning. R and Z are the r, z coordinates of the point of interest, RZ is the radius of the disc, or of the cross section of the cylinder, RL is the length of the cylinder for the finite cylinder case and RS is the output from this subroutine and represents the following dimensionless form of the integral in equation (3), fo RS =
Yl(t)Jdxt)e-" dt
(4)
where x = r/a and y = zla. In all practical applications, one would be interested in calculating gz at a large number of points for which experimental observations are planned or already available. Hence a short main program somewhat along the lines of the test program which we present in the Appendix would have to be written to define the points of interest. This main program must have the common declaration (~OMMON X, I/, ISHAPE and ISHAPE must be defined before reference is made
315
to the subroutine AXIAL. This parameter ISHAPE must be assigned a value of 1, 2, or 3 depending on whether calculations are to be made for a circular disc, an infinite cylinder, or a finite cylinder respectively. The output RS returned by the subroutine AXIAL may be multiplied by the appropriate units in the main program to give g~ at the required point. The subroutine AXIAL calls two functions PINT and FINK. The function FINK performs the task of calcularing the integrand in equation (4) for arbitrary values of T for a given X and Y. The X and Y are referenced through the blank common mentioned previously. The function RINT to which a call is made by the statement RS = RINT (FINK, A, B, EPS) integrates the function from A to B by an adaptive Simpson's quadrature. EPS is an error parameter, and can be set to 0.0001. The lower limit for the integral A is always set to 0, and B to 15.0/¥, for most situations. For very small values of Y, the integration is bypassed and the result is set to the analytically known result for Y = 0. tmULTS~OM nOGRAM The main aim of this article is to provide a program which calculates gz at the point (r, z) for objects with radial symmetry, thus reducing the labor of looking into tabulated functions and the interpolations. We therefore have used the program to obtain illustrative results for special situations, for which exact results have been reported earlier in the literature. As a first application, we have calculated the gz at various (r, z) for a circular disc. Our test program results reproduce exactly the results as presented by Singh (1977a, table 1). We did calculations at a number of other points for display from GRPLT with a line printer as well as with the CALCOMP plotter. Our results are somewhat at variance with those reported by Parasnis (Parasnis, 1961), at times by about 5 percent. If anything, this implies that his truncation at the fourth term yields reasonable results. Our direct integration technique is more versatile, convenient and accurate than an expansion in terms of Legendre functions. We also have used the program for the situation of an infinite cylinder, and finite cylinder as well as for an hemisphere and a cone. For points on the axis, the results agree closely with those from the analytic ones for these situations. For a compilation of the appropriate expressions the reader is referred to the review article by Talwani (1973). ADArrA~oNsANDEXTENSIONS The program presented here is applicable as such for a large number of situations of practical intero, t. With a suitably written program segment it may be used for any composite shape which can be built up of cylinders, cones, and spheres. By allowing for different densities of the discs into which the three-dimensional body may be divided, the program may be extended to take account of the axial variations in the density.
316
M.M. PANTand J. G O v m D ~
~,ud£NCES AMamowitz, M., and Stegun, I. A., 1964, Handbook of mathcmatical functions: Dover Publ., Inc., New York, p. 369-370. Bateman, H., 1959, Partial differential equations of mathematical physics: Cambridge Univ. Press, Cambridge, 469p. Milsom, J., and Worthington, G. A., 1977, Computer pmgralns for rapid computation of gravity effects of two-dimensional and three-dimensional bodies: Computers & Geosciences, v. 3, no. 2, p. 269-281. Nabiighian, M. N., 1962, The gravitational attraction of a right vertical circular cyfinder at points external to it: Geof. pura e Appl. v. 53, p. 45-51. Nettleton, L. L., 1940, Geophysical prospecting for oil: McGraw Hill Book Co., New York, 444 p. Nettleton, L. L., 1942, Gravity and magnetic calculations: Oeophysics, v. 7, no. 3, p. 293--310. Parasnis, D. S., 1961, Exact expressions for the gravitational
attraction of a circular lamina at all points of space and of a fight circular vertical cylinder at points external to it: Oeophys. Prosp., v. 9, p. 382-398. Rao, B. S. R., and Radhakrishnamurty, 1966, An expression for the gravitational attraction of circular plate: Indian Jour. Pure Appi. PhD., v. 4, p. 276--278. Singh, S. K., 1977, Gravitational attraction of a circular disc: Geophysics, v..42, no. 1, p. 111-113. Singh, S. K., 1977a, Gravitational attraction of a vertical right circular cylinder! Geophys. Jour. Roy. Astr. Sc~., v. 50, p. 243-246. Talwani, M. 1973, Computer usages in the computation of gravity anomalies,/n Methods of computational physics, v. 13: ed. by B.A. Bolt, Academic Press, New York, p. 343--389. Wolfe, P., 1971, Eigentunctions of the integral equation for the potential of the charged disk: Jour. Math. Phys., v. 12, p. 1215-1222.
AEPENI)C(
PROGRAM TEST(INPUT,OUTPUT,TAPES=INPUT,TAPE6=OUTPUT,PUNCH) COMMON X,Y, ISHAPE DIMENSION S X C l l 0 ) , S Y C l l 0 ) , S F B C I 1 0 , 1 1 0 ) , S P H C I 0 0 ) , D V ( 5 0 , ~ ) C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
THIS IS A TEST PROGRAM WHICH IS DESIGNED TO SHOW HOW THE SUBROUTINE AXIAL AND THE RELATED ROUTINES ARE TO BE ACCESSED IN ORDER TO CALCULATE THE Z COMPONENT OF THE GRAVITATIONAL F I E L D . THERE ARE FOUR SECTIONS OF THIS PROGRAM WHICH DEALS WITH THE CALCULATIONS FOR OBdECTS SUCH AS A D I S C , I N F I N I T E CYLINDER,HEMISPHERE AND A RIGHT CIRCULAR CONE RESPECTIVELY. THE VARIABLE ISHAPE IN THE BLANK COMMON IS TO BE SET BY THE USER IN THE MAIN PROGRAM ACCORDING TO THE FOLLOWING DESCRIPTION. ISMAPE
SX SY NX NY
= 1 2 3
FOR CIRCULAR DISC I N F I N I T E CYLINDER F I N I T E CYLINDER
ARRAY OF R COORDINATES ARRAY OF Z COORDINATES NUMBER OF ELEMENTS IN THE SX ARRAY NUMBER OF ELEMENTS IN THE SY ARRAY
NDISC NUMBER OF DISCS INTO WHICH THE OBdECT IS TO BE D I V I D E D . NOTE THAT THIS IS NOT NECESSARY IN THE CASE OF A DISC OR A CYLINDER. IN THE TEST EXAMPLE DESCRIBED HERE NDISC HAS BEEN ARBITRARILY SET TO l O 0 . DSX AND DSY DE4OTE THE INTERVAL BETWEEN TWO SUCCESSIVE VALUES OF SX AND SY RESPECTIVELY. IPLT .GT. .LE.
CONTROLS THE PLOT ON THE L I N E 0 0
PRINTER
PLOT IS PRODUCED PLOT 1S OMITTED
FOR THE CASE OF A RIGHT CIRCULAR CONE THE FOLLOWING THREE PARAMETERS MUST BE GIVEN AS INPUT. A BETA ZER
RADIUS OF THE BASE OF THE CONE HALF THE ANGLE OF THE CONE I N DEGREES DENOTES Z0 IN THE FIGURE AND CORRESPONDS TO THE HEIGHT OF THE APEX OF THE CONE ABOVE THE Z=0 PLANE
ITEST CONTROL VARIABLE FOR THIS TO PERFORM S P E C I F I C TASKS =l 2
C
3
C C C C C
~
TEST PROGRAM ONLY
DISC I N F I N I T E CYLINDER HEMISPHERE CONE
FOR THE CASE OF A HEMISPHERE, THE RADIUS OF THE HEMISPHERE IS TO BE GIVEN AS THE INPUT. ( REPRESENTED BY THE VARIABLE RHS )
Compu~tionofgmvitye~ectsoftwo-~d ~e~imension~bodies
317
C C C C C C C C
THE OUTPUT FROM THE SUBROUTINE AXIAL IS CONTAINED IN THE VARAIBLE RS. THE VARIABLE SFB IS A MATRIX WHICH IS USED TO COPY THE RESULT AS A FUNCTION OF RZERO AND Z RESPECTIVELY. WRITE(6,10) 10 FORMAT(IHI) 20 CONTINUE READ(5,30) ITEST,NX,NY,IPLT,ND1SCpDSX,DSY 30 F O R M A T ( S I S , 2 F I 0 . 5 )
C C C C C
THE FOLLOWING STATEMENT CHECKS FOR THE END OF F I L E ON INPUT AND IS APPLICABLE TO ONLY CONTROL DATA MACHINES. FOR OTHER COMPUTERS CHECK THE FORTRAN LANGUAGE MANUAL.
IF(EOF(5)) 3 1 0 , ~ 0 W0 CONTINUE WRITE(6,SD) ITEST,NX,NY,IPLT,NDISC,DSX,DSY 50 FORMAT(IHO,SHITEST,IX,I),2X,3HNX ,I5,1X,)HNY , I 5 , 1 X , S H 1 P L T x 6HNDISC ,15,1X,3flDSX,FI0.5,1X,3HDSY,FI0.5) C C C
GENERATE
60 70
80
90 100
C C C C C
,15,1X,
SX AND SY ARRAYS
IF(NDISC.LE.0) NDISC =100 SX(1)=0.0 SY(1)=0.0 DO 60 I=2,NX SX(1)=SX(I-I)+OSX CONTINUE DO 70 I =2,NY SY(1)=SY(I-I)+DSY CONTINUE GO TO (80, I~0, 160, 210) ITEST CONTINUE ISHAPE = 1 DO 100 I = 1,NX DO 90 d = 1,NY R = SX(1) Z = SY(J) CALL A X I A L ( R , Z t l . 0 , 1 . 0 , R S ) SFB(I,J) = RS CONTINUE CONTINUE IF(.NOT.(IPLT.GT.0)) GO TO 280 WRITE(6,10) M = 50 OUTPUT PLOT IS WRITTEN ON UNIT 6 . IF THE PLOT IS DESIRED ON ANY OTHER UNIT , CHANGE NC TO THE APPROPRIATE VALUE
NC:6 N=k DO 130 d = l , N DO 120 I = I , M D V ( l , d ) = SFB(I,J)ISFB(I,d) 120 CONTINUE 130 CONTINUE CALL G R P L T ( N C , S X , D V ~ M p N ) GO TO 280 1W0 CONTINUE C C C
INFINITE
CYLINDER
DO 155 I = I , N X OO 150 J = I , N Y R = SY(1) Z = SY(d) ISHAPE = 2 CALL A X I A L ( R , Z , I . 0 , 1 . 0 , R S ) SFB(I,J)=RS 150 CONTINUE 1SS CONTINUE GO TO 280 160 CONTINUE C
¢
318
M.M. PANTand J. GOVINDAP~JAN C C C
HEMISPHERE READ THE RADIUS OF THE HEMISPHERE 165
170
180 190 200 210 C C C
READ(5,165) RHS FORMAT(F10.5) NDM = NDISC - 1 DO 200 KK = I,NX R = SX(KK) DO 190 K = 1,NY XF = SY(K) DX=RHS/FLOAT(NDISC) DO 170 I = I , N D I S C XX=(l-l)xOX RZ=SQRTCXXX(2.OXRHS-XX)) Z=XF+RHS-XX ISHAPE = 1 CALL A X I A L ( R , Z , R Z , I . 0 , R S ) SPH(1) = RS CONTINUE ZPP = 0 . 5 ; : ( S P H ( 1 ) + SPH(NDISC))XDX OO 180 I=IpNOM ZPP=ZPP+SPHCI~XDX CONTINUE SFB(KK,K) = ZPP CONTINUE CONTINUE GO TO 280 CONTINUE RIGHT CIRCULAR CONE
READ(S,220) A, BETA,ZER 2"20 FORMAT(3F10.5) WR1TEC6,230) A,BETA,ZER 230 FORMAT(IH0,1OX,3H A = , F 6 . 2 , 2 X , 6 H C C C
BETA=,2X,FI0.2,2X,SH ZER=,FL0.3,/)
CONVERTBETA FROM DEGREES TO RADIANS PI=q.0XATAN(I,0) BETA=BETA;PI/180.0 NDM = NDISC - 1 CT=I.0/TANCBETA) DX=(AXCT-ZER)/NOISC DO 270 1 = I,NX DO 260 J = 1,NY R : SX(l) DO'2~0 K = 1,NDISC RZ =( ZER + K~DX)XTAN(BETA) Z = SY(J) + KXOX ISHAPE = 1 CALL A X I A L ( R , Z , R Z , I . 0 , R S ) 2~0 SPH(K)=RS SUM=0.5~DXX(SPH(1)+SPH(NDISC)) DO 250 M=2,NDM SUM=SUM+OXXSpH(M)
250 CONTINUE S F B ( 1 , J ) = SUM 260 CONTINUE 270 CONTINUE 280 CONTINUE WRITE(6,10) DO 285 I=I,-NX WRITE(6,290) C S F B C I , J ) , J = I , N Y ) 285 CONTINUE 290 FORMATCIHO,CSF15.6,/)) PUNCH 3 0 0 , ( ( S F B ( 1 , J ) , I = I , N X ) , J = I , N Y ) 300 FORMAT(8F10.6) GO TO 20 310 CONTINUE STOP END SUBROUTINE A X I A L ( R , Z , R Z , R L , R E S ) COMMON X , Y , ISHAPE EXTERNAL FINK C C C C
THIS SUBROUTINE IS THE KEY SUBROUTINE WHICH ACTS AS THE INTERFACE BETWEEN THE USER DEFINED MAIN PROGRAM AND THOSE ROUTINES WHICH PERFORM THE ACTUAL CALCULATION
Computation of gravity effects of two- and three-dimensionalbodies C C C C C C C C C C C C C C
319
INPUT . . . . . . . . . . . R AND Z ARE RESPECTIVELY THE COORDINATES OF THE POINT AT WHICH THE F I E L D IS TO BE CALCULATED RZ
THE RADIUS OF THE DISC OR THE CYLINDER
RL THE LENGTH OF THE F I N I T E
CYLINDER
OUTPUT . . . . . . . . . RES
VARIABLE
WHICH CONTAINS THE RESULT
YZ = 0 . 0 0 1 X=R/RZ Y=Z/RZ A=0.0 B = 20.0/Y EPS=0.0001 IF(Y.LT.YZ) GO TO 10 IF (Y . L T . 0 . 0 1 5 ) EPS = 0 . 0 0 0 0 1 RES=Ri N T ( F I N K , A, B, EPS) IF(ISHAPE.LE.2) GO TO 20
Y=(Z+RL)/RZ A = 0.0 B = 20.0/Y EPS = 0 . 0 0 0 1 RESA=RI NT(F I N K , A , B , E P S ) RES=RES - RESA GO TO 20 10 CONTINUE IF(X.LT.1.0) RES = 1 . 0 IF(X.EQ. 1.0) RES=0.5 IF(X.GT.1.0) RES=0.0 20 CONTINUE RETURN END
FUNCTION R I N T ( F , A , B t E P S ) DIMENSION D X C 3 0 ) , E P S P C 3 0 ) , X 2 C 3 0 ) , X 3 C 3 0 ) , F 2 C 3 0 ) , F 3 ( 3 0 ) , F ~ ( 3 0 ) ,
=
FMP(~0),FBP(~0),EST2(~O),EST3(~O),PVAL(~0,3) INTEGER RTRNC30)
C
C C C
THIS FUNCTION SUBPROGRAM PERFORMS THE INTEGRATION BY A NONRECURSIVE ADAPTIVE SINP$ONS INTEGRATION FORNULA. THE FORTRAN CODE IS BASED ON THE ALGORITHM 182
C C
PUBLISHED
IN THE COMMUNICATIONS OF ACM, PAGE 2~Np
196k.
LVL=0 ABSA=I.0 EST=I.0 DA=B-A FA=F(A) FM=~.0xFCCA+B)/2.0) FB=F(B) 10 CONTINUE LVL=LVL+I DX(LVL)=OA/3.0 SX=DX(LVL)/6.0 FI=~.0XF(A+DX(LVL)/2.0) X2(LVL)=A+OX(LVL) F2(LVL)=F(X2(LVL)) X3CLVL)=X2CLVL)+DXCLV L) F3(LVL)=F(X3(LVL)) EPSP(LVL)=EPS Fk(LVL)=h.0;:F(X3(LVL)+OX(LVL)/2.0) FMP(LVL)=FM ESTI=(FA+FI+F2(LVL));SX FBP~LVL)=FB EST2(LVL):(F2(LVL)+F3(LVL)+FM)XSX EST3(LVL)=(F3(LVL)+F~(LVL)+FB)XSX SUM=ESTI+EST2(LVL)+EST3(LVL) ABSA=ABSA-ABS(EST)+ABS(EST1)+ABS(EST2(LVL))+ABS(EST3(LVL))
IF(.NOT.((ABS(EST-SUM).LE.EPSP(LVL)XABSA).AND.(EST.NE.1.0)) X(LVL.GE.30))GO 20 CONTINUE LVL=LVL-I
TO 30
.OR.
3~
M.M. PAN'r~dJ.OOV[NIDA~JAN
30
~0
50
60
PVAL(LVL,RTRN(LVL))=SUM GO TO ( ~ 0 , 5 0 , 6 0 ) RTRN(LVL) CONTINUE RTRNCLVL)=I DA=DX(LVL) FM=F1 FB=F2(LVL) EPS=EPSPCLVL)/1.7 EST=ESTI GO TO 10 CONTINUE RTRNCLVL)=2 DA=DX(LVL) FA=F2(LVL) FM=FHP(LVL) FB=F3(LVL) EPS=EPSP(LVL)/I.7 EST=EST2(LVL) A=X2CLVL) GO TO 10 CONTINUE RTRNCLVL)=3 DA=DX(LVL) FA=F3(LVL) FM=F~(LVL) FB=FBP(LVL} EPS=EPSP(LVL) EST=EST3CLVL) A=X3(LVL) GO TO 10 CONTINUE SUM=PVAL(LVLtl)+PVAL(LVL,2)+PVAL(LVL,3) | F C L V L . G T . I ) GO TO 20 RINT=SUM RETURN END FUNCTION FINKCT)
C C C C C C
THIS FUNCTION SUBPROGRAM CALCULATES THE BESSEL FUNCTIONS BASED ON THE FORMULAE AS GIVEN IN THE HANDBOOK OF MATHEMATICAL FUNCTIONS BY ABRAMOWITZ AND STEGUN. DIMENSION A ( 6 ) , B ( 7 ) , C ( 7 ) , D ( 6 ) , E ( 7 ) , F ( 7 ) COMMON X , Y , ISHAPE DATA A I - 2 . 2 ~ 9 9 9 9 7 , 1 . 2 6 5 6 2 0 8 , - 0 . 3 1 6 3 8 6 6 , 0 . 0 ~ k ~ 7 9 , - 0 . 0 0 3 9 ~ , 1 0.0002100/
C DATA B / 0 . 7 9 7 8 8 ~ 5 6 , - 0 . 0 0 0 0 0 0 7 7 , - 0 . 0 0 5 5 2 7 h 0 , - 0 . 0 0 0 0 9 5 1 2 ,
1
0.0013727,-0.00072805,0.0001~k761
C 1
DATA C / - 0 . 7 8 5 3 9 8 1 6 t - 0 . 0 ~ 1 6 6 3 9 7 , - 0 . 0 0 0 0 3 9 5 ~ , 0 . 0 0 2 6 2 5 7 3 , -0.0005~125,-0.000293~3,0.00013558/
1
DATA D / - 0 . 5 6 2 k 9 9 8 5 , 0 . 2 1 0 9 3 5 7 ~ , - 0 . 0 3 9 5 ~ 2 8 9 , 0 . 0 0 ; ~ 3 3 1 9 , - 0 . 0 0 0 3 1 7 6 1 , 0.00001109/
1
DATA E 10.79788~56,0.00000156,0.01659667,0.00017105,-0.002~9511, 0.00113656,-0.00020033/
1
DATA F / - 2 . 3 5 6 1 9 k k 9 , 0 . 1 2 k 9 9 6 1 2 , 0 . 0 0 0 0 5 6 5 0 , - 0 . 0 0 6 3 7 8 7 9 , 0 . 0 0 0 7 ~ 3 h B , 0.0007982~,°0.00029166/
C C C C DATA TL I I . O E - 0 7 / C C C
MAIN LOOP GENERATING BESSEL FUNCTIONS I F C T . L E . T L ) T = TL XT=X~T yT=Y~:T IF(T-3.0) 10,10,30 10 CONTINUE X3=T/3.0 X32=X3~X3 FOI=0.5 DO 20 g = l , 6 Fdl=FJ1+D(K)=X32~;K 20 CONTINUE
Computation of gravity effects of two- and three-dimensional bodies
30
~0 50 60
70 80
90 100
Fdl=FdlXT GO TO 50 CONTINUE X3=3.0/T FFI=E(1) FF2=F(I) DO ~0 K = 2 , 7 X3K=X3X=C(K-1) FFI=FFI+E(K)XX3K FF2=FF2+FCK)XX3K CONTINUE FdI=FFI=tCOS(FF2+T)/SQRT(T) CONTINUE IF(XT-3.0) 60,60,80 CONTINUE X3=XT=XT/9.0 FdZ=I.0 DO 70 K = l , 6 FJZ=FdZ+A(K)=X3=XK CONTINUE GO TO 100 CONTINUE X3=3.0/XT FFI=B(I) FF2=C(1)+XT DO 90 K = 2 , 7 X3K=X3X~(K-1) FFI=FFI+B(K)=X3K FF2=FF2+C(K)~X3K CONTINUE FdZ=FFI=COSCFF2)/SQRT(XT) CONTINUE FINK=FdZ=Fdl~EXP(-YT) IF ( I S H A P E . E Q . 2 ) F I N K = F I N K / T IF ( I S H A P E . E Q . 3 ) F I N K = F I N K / T RETURN END SUBROUTINE GRPLT(NC,X,YsM,N)
C C C C
THIS SUBROUTINE PERFORMS THE TASK OF CREATING A PLOT ON THE LINE PRINTER.
C
NC IS THE NUHBER OF THE DATA CHANNEL
C C C C C
THI$ SUBROUTINE IS BASED ON THE ALGORITHM PUBLISHED IN THE COMMUNICATIONS OF ACM, P ~ 9 2 , 1 9 7 1 AS ALGORITHM ~ 1 2 .
C C C
USUALLY UNIT 6
DIMENSION X(M),Y(H,N),LP(50),IW(IOO),PS(IO0) INTEGER S,LsPS DATA X E R R / 0 . 0 0 0 I / , Y E R R / 0 . 0 0 0 1 / DATA L / 5 0 / , 5 / 8 5 / DATA I B / 1 H / DATA MC/1H+/ DATA I C 1 / I H A / DATA I C 2 / I H B / DATA I C 3 / 1 H C / DATA IC~/1HO/ FIND XHAX AND XMIN XMAX=X(I) XMIN=X(1) DO 10 l=2,N IF(X(1).GT.XMAX) IF(X(I).LT.XMIN) 10 CONTINUE
C C C
-
XMAX=XCI) XMIN=X(I)
FIND YHAX AND YHIN
YMIN=YCI,I) YMAX=Y(1,1) DO 50 DO 20
I=1,~ d=l,N
IF(Y(I,d).GT.YMAX) YMAX=YCI,d) IF(Y(I,d).LT.YMIN) YMIN=Y(I,J) 20 CONTINUE
321
M.M. PANT and J. GOVlND~AN
322
30 CONTINUE 00 ~0 I=I,M OO 35 d = l , N Y(I,d)=T(I,d)/YMAX 35 CONTINUE 40 CONTINUE XMIN=0.0 YMIN=0.0 YMAX=I.O IFCCABS(XMAX-XMIN).LT.XERR).OR.CABS(YMAX-YMIN).LT.YERR))RETURN
PX=FLOAT(L-1)/(XMAX-XMIN)
50 60
70 80 85
PY=FLOAT($-I)I(YMAX-YMIN) DO 50 I = I , L LP(I)=0 CONTINUE KL=M+I CONTINUE KL=KL-1 IF(KL.LT.I) GO TO 70 KR=I+INT(CX(KL)-XMIN)XPX+0.5) LP(KR)=KL GO TO 60 CONTINUE WRITE(NC,80) FORMAT(1HO) WRITECNC,$5) FORMATC1H , ~ O X , I S H NORMALIZED GZ ) DO 90 d = l , S
IW(J)=MC 90 CONTINUE WRITE(NCt95) 95 F O R M A T ( l O X , 3 H O . 0 , 8 2 X , ) H I . 0 ) WRITECNC,100) ( I W ( d ) , d = I , S ) 100 FORMAT(1H , 1 0 X , 1 0 0 A 1 ) DO 250 I = I , L PS(1)=l PS(S)=I MS=S-1 DO 110 J=2,M$ PS(d)=2 110 CONTINUE IF(.NOT.(LP(I).GT.0)) GO TO 150 120 CONTINUE IF(.NOT.(N.EQ.1)) GO TO 130
INTR=I+INT(0.5+PY::(Y(LP(1),I)-YMIN)) PS(INTR)=3 GO TO 150 130 CONTINUE d=N+l 140 CONTINUE
J=d-I IF(J.LT.I) GO TO 150 INTR=I*INT(0.5+PY;:(Y(LP(1),J)-YMIN)) PSCINTR)=J+2 GO TO 14O 150 CONTINUE DO 220 J = l , S
IW(U)=IB IF(.NOT.(PS(d).EQ.1))
GO TO 170
IW(J)=MC GO TO 220 170 CONTINUE IF(.NOT.(PS(J).EQ.2))
GO TO 180
lW(d)=I8 GO TO 220 180 CONTINUE IFC.NOT.(PSCd).EQ.3))
GO TO 190
IW(d)=ICl GO TO 220 190 CONTINUE IF(.NOT.(PS(J).EQ.~))
GO TO 200
GO TO 220 200 CONTINUE IF(.NOT.(PS(d).EQ.5))
GO TO 210
IW(J)=IC2
IW(d)=IC3 GO TO 220 210 CONTINUE
IW(d)=IC~
Computation of gravity e~ects of two- and three-dimensional bodies
220 CONTINUE IFC.NOT.CI.EQ.1)) GO TO 230 WRITECNC,225) (IWCd),d=l,S) 225 FORMATC6Xt3H0.0t2X,100A1) GO TO 250 230 CONTINUE IFC.NOT.CI.EQ.25)) GO TO 2k0 WRITECNC,235) (IWCd),J=I,S) 255 FORMAT(6X,3HR/A,2X,10OA1) GO TO 250 2k0 CONTINUE WRITECNC,2qS) ( I W ( d ) e d = l , $ ) 2k5 FORIqATC1H ,10Xpl00Al) 250 CONTINUE 260 FORMAT(1H ,11) DO 270 d=l,S IW(d)=MC 270 CONTINUE WRITECNCe275) XMAX,CIWCJ),J=lt$) 275 FORRAT(6X,Fk.2,1X,10OAI) WRITE(NC,260) WRITECNC;2$0) XMINtXMAX,YMIN,YMAX 280 FORHAT(1HOpIOX,kHXMIN,FIO.k,11X,~HXMAX,F10.~,llX,kHYMIN,FlO.qp 111X,~HYMAJ(,FI0.N) RETURN END
323