DECLARE SUB pict1 () DECLARE SUB Pict2 () DECLARE SUB Pict3 () DECLARE SUB Pict4 () DECLARE SUB Square (a!, b!, c!, d!) DECLARE SUB Pager () REM GE Copeland REM Physics Department REM Old Dominion University REM Norfolk, VA 23529 REM Udated 3/2001 86 PRINT CHR$(12); PRINT "Go directly to the calculation? Y=1" INPUT Jump IF Jump = 1 THEN 11102 88 REM 100 REM for further information about this program contact 200 REM conduit p.o. box 388 iowa city, iowa 52240 (319)353-5789 300 REM registry number: phy003 400 REM last revision 5/77 by neil s. ferguson for conduit 500 REM this program is distributed with the support of the national 600 REM science foundation grant no. Sed75-06596. Any opinion,findings, 700 REM conclusions, or recommendations expressed or implied are those 800 REM of the authors and do not necessarily reflect the views of the 900 REM national science foundation 1000 REM copyright by houghton-mifflin, 1976. All rights reserved. 1100 REM "Poisso" -- (basic program begins at line 290) 1200 REM 1300 REM 1400 REM description: 1500 REM 1600 REM this program solves poisson's equation using simplest 1700 REM central difference equation. 1800 REM 1900 REM 2000 REM 2100 REM instructions: 2200 REM 2300 REM this program is taken from "using computers in 2400 REM physics" by john r. merrill and should be used in conjunction 2500 REM with that text. The programs are distributed by project 2600 REM conduit. CALL pict1 CALL Pager CALL Pict2 CALL Pager CALL Pict3 CALL Pager CALL Pict4 11102 REM conv = .0005 11300 PRINT "Have a good time - output will be a surface plotting file in the" PRINT "proper format for WPLOT. Convergence Criterion is for the fractional" PRINT "change over an iterations to be less than"; conv; "." 11400 PRINT 11500 PRINT "Do you wish to choose the potentials v1,v2,v3 and v4" 11600 PRINT "(Yes or No) - No will automatically supply some values."; 11700 INPUT p$ 11800 IF LEFT$(p$, 1) = "n" OR LEFT$(p$, 1) = "N" THEN GOTO 13900 11900 GOTO 12300 12000 REM shouldn't be here 12100 PRINT "I don't understand your selection." 12200 GOTO 11400 12300 PRINT 12400 PRINT "What values do you want for v1,v2,v3 and v4 (all less than 30 volts)"; 12500 INPUT a, b, c, d 12600 IF a > 30 THEN GOTO 13500 12700 IF a < 0 THEN GOTO 13500 12800 IF b > 30 THEN GOTO 13500 12900 IF b < 0 THEN GOTO 13500 13000 IF c > 30 THEN GOTO 13500 13100 IF c < 0 THEN GOTO 13500 13200 IF d > 30 THEN GOTO 13500 13300 IF d < 0 THEN GOTO 13500 13400 GOTO 15110 13500 PRINT 13600 PRINT "This program is designed to be relatively rapid and simple;" 13700 PRINT "thus your selections must be between 0 and 30." 13800 GOTO 12300 13900 PRINT 14000 PRINT "I will let v1=1, v2=10, v3=19 and v4=10." CALL Pager 14100 a = 1 14200 b = 10 14300 c = 19 14400 d = 10 CALL Square(a, b, c, d) CALL Pager 15110 REM program starts 15510 REM array size i1 = 10 j1 = 10 DIM V(i1, j1), u(i1, j1) 15300 v1 = 0 15400 n = 0 15500 FOR i = 1 TO i1 FOR j = 1 TO j1 V(i, j) = 0! NEXT j NEXT i 15800 REM setting v(i,j) boundary conditions 15900 V(1, 1) = a 16000 V(i1, 1) = b 16100 V(i1, j1) = c 16200 V(1, j1) = d 16300 REM slopes of straight lines connecting the corners 16400 m1 = (b - a) / (i1 - 1) 16500 m2 = (c - b) / (j1 - 1) 16600 m3 = (c - d) / (i1 - 1) 16700 m4 = (d - a) / (j1 - 1) 16800 REM calculate the potentials on the four sides 16900 REM side ab and side dc 17000 FOR i = 2 TO i1 - 1 17100 V(i, 1) = V(1, 1) + m1 * (i - 1) 17200 V(i, j1) = V(1, j1) + m3 * (i - 1) 17300 v1 = v1 + V(i, 1) + V(i, j1) 17400 NEXT i 17500 REM side ad and bc 17600 FOR j = 2 TO j1 - 1 17700 V(1, j) = V(1, 1) + m4 * (j - 1) 17800 V(i1, j) = V(i1, 1) + m2 * (j - 1) 17900 v1 = v1 + V(1, j) + V(i1, j) 18000 NEXT j 18100 v1 = v1 / (2! * (i1 + j1)) 18200 FOR i = 2 TO i1 - 1 18300 FOR j = 2 TO j1 - 1 18400 V(i, j) = v1 18500 NEXT j 18600 NEXT i 18700 v2 = 0 18800 FOR i = 1 TO 10 18900 FOR j = 1 TO 10 19000 u(i, j) = V(i, j) 19100 NEXT j 19200 NEXT i 19300 FOR i = 2 TO i1 - 1 19400 FOR j = 2 TO j1 - 1 19500 r = 1 19600 V(i, j) = (u(i + 1, j) + u(i - 1, j) + u(i, j + 1) + u(i, j - 1) + r) / 4 19700 IF v2 > ABS((V(i, j) - u(i, j)) / V(i, j)) THEN GOTO 19900 19800 v2 = ABS((V(i, j) - u(i, j)) / V(i, j)) 19900 NEXT j 20000 NEXT i 20100 n = n + 1 20200 IF v2 > conv THEN GOTO 18700 20210 PRINT 20300 PRINT "# iterations over grid ="; n 20310 PRINT 20312 PRINT "The data is read in the following manner:" 20314 PRINT " v(1,10) ,v(2,10),. . . . . . . . . . .v(10,10)" 20320 PRINT " v(1,9) ,v(2,9). . . . . . . . . . . .v(10,9)" 20322 PRINT " v(6,9), . . . . " 20324 PRINT 20326 PRINT " v(1,1)=, . . . . . . . . . . . . . . .v(10,1)" PRINT 20400 FOR j = j1 TO 1 STEP -1 20500 FOR i = 1 TO i1 20600 PRINT USING "##.#### "; V(i, j); 20700 NEXT i 20800 PRINT 20900 NEXT j REM find max and min V Vmax = -9.999999E+27 Vmin = 1E+29 FOR j = 1 TO j1 FOR i = 1 TO i1 IF Vmax < V(i, j) THEN Vmax = V(i, j) IF Vmin > V(i, j) THEN Vmin = V(i, j) NEXT i NEXT j PRINT "The max value of V in Volts is "; Vmax; "and the min value is "; Vmin CALL Pager PRINT " The PLOT.PLT file is being created, return to windows and run WPLOT." PRINT "You will see three 3D surface plots." REM now write out the plotting file OPEN "PLOT.PLT" FOR OUTPUT AS #2 REM example PRINT #2, "CLRDATA" PRINT #2, "CLRPLOT"; "" PRINT #2, "CLOSE" PRINT #2, PRINT #2, "TITLE 3-D Plot sliced "; j1; "x"; i1; " grid" PRINT #2, "ZLABEL Y" PRINT #2, "XLABEL X" PRINT #2, "YLABEL Potential in Volts " PRINT #2, "YSTART "; Vmin PRINT #2, "YEND "; Vmax PRINT #2, "YSTEP "; INT((Vmax - Vmin) / 20!) PRINT #2, "XSTART 1" PRINT #2, "XEND "; i1 PRINT #2, "XSTEP 1" PRINT #2, "XTICS 3" PRINT #2, "ZSTART 1" PRINT #2, "ZEND "; j1 PRINT #2, "ZSTEP 1" PRINT #2, "ZTICS 0" PRINT #2, "READZ "; i1 FOR k1 = 1 TO i1 PRINT #2, k1; NEXT k1 PRINT #2, PRINT #2, "READTAB "; i1; " "; i1 FOR j = 1 TO j1 STEP 1 PRINT #2, j; FOR i = 1 TO i1 PRINT #2, USING "##.#### "; V(i, j); NEXT i PRINT #2, NEXT j PRINT #2, "PLOT3D" PRINT #2, PRINT #2, "TITLE Rotated 3-Dimensional Plot" PRINT #2, "ROTATE" PRINT #2, "PLOT3D" PRINT #2, PRINT #2, "TITLE Surface Plot from Poisson" PRINT #2, "PLOT3D SURFACE" CLOSE #2 21000 END SUB Pager PRINT , "Push return to continue"; INPUT l$ PRINT CHR$(12) END SUB SUB pict1 2700 PRINT "Title: POISSON--a program for finding a solution to Poisson's" 2800 PRINT "equation." 2900 PRINT 3000 PRINT 3100 PRINT "This program is designed to demonstrate a computer method" 3200 PRINT "for solving Poisson's equation. It uses a technique similar" 3300 PRINT "to that used in program LAPLACE.BAS (the student should run" 3400 PRINT "LAPLACE.BAS for background information concerning Laplace's" 3500 PRINT "and Poisson's equations and the numerical method used here" 3600 PRINT "and there). You use the same ideals to approximate del**2 V," PRINT "where del is the Nabla Operator. In Cartesian coordinates, it is" PRINT PRINT " 2 2 2" PRINT " { d + d + d } V(x,y,z) = -p(x,y,z)" PRINT " ----- ------ ----- --------" PRINT " dx^2 dy^2 d z^2 Eo" PRINT PRINT " and Eo is the permitivitty of free space, and P is the charge" PRINT "density ." END SUB SUB Pict2 3800 PRINT "The difference equation becomes (run laplac.bas if you" 3900 PRINT "don't understand some of the notation used):" 4000 PRINT 4100 PRINT "V(x,y)=(approx.) (V(x+d,y)+V(x,y+d)+V(x-d,y)+V(x,y-d)+" 4200 PRINT " (D**2/a)p(x,y))/4" 4300 PRINT , "where Eo=epsilon and p=charge density - in 2D!" 4500 PRINT "The program Poisson illustrates the technique in a rec-" 4600 PRINT "tangular region. For simplicity p/a is chosen to be 1" 4700 PRINT "everywhere. Hence Poisson's equation is of the form" 4800 PRINT "del^2 V(x,y,z)=-p/a=-1 (charges are present)." END SUB SUB Pict3 7500 PRINT "The solution to Laplace's and Poisson's equations in" 7600 PRINT "multiply connected regions (regions with holes) is more" 7700 PRINT "complicated. But in some cases the problem can be solved" 7800 PRINT "(additional information can be found in the library under" 7900 PRINT "numerical methods--or a particularly good reference is:" 8000 PRINT "Applied Numerical Methods by Luther and Wilkes). Such" 8100 PRINT "multiply connected region problems are common." PRINT "An example" 8200 PRINT "is:" 8300 PRINT 8400 PRINT "Place a point charge in an otherwise charge-free region which " 8500 PRINT "has fixed potentials on its boundries. One method of so-" 8600 PRINT "lution is to fix the potential on the grid points near" 8700 PRINT "the point charge as being due entirely to the point charge." END SUB SUB Pict4 8800 PRINT "You then solve Laplace's equation for a region having a hole " 8900 PRINT "in it--the hole is the region around the point charge" 9000 PRINT "you must still assume that all the values of V on the" 9100 PRINT "boundry are fixed, even those right around the hole." PRINT 9300 PRINT "Now in this particular program you will be given a square" 9400 PRINT "region of space - 10x10" 9500 PRINT 9600 PRINT " v4=?-----------v3=?" 9700 PRINT " 1 1" 9800 PRINT " 1 1" 9900 PRINT " 1 1" 10000 PRINT " 1 1" 10100 PRINT " 1 1" 10200 PRINT " v1=?-----------v2=?" 10300 PRINT 10400 PRINT "with uniform charge density p=a throughout the region" 10410 CALL Pager 10900 PRINT "Our Grid will be a 10x10 matrix and will compute it over 100" PRINT "times..... takes time to do this!" 11000 PRINT "At the start of the program you will be given the option" 11100 PRINT "of choosing v1(1,1), v2(10,1), v3(10,10) and v4(1,10)." 11200 PRINT " These are the potentilas at the four corners of the square." PRINT " and there is a constant space charge on the 2D-surface." END SUB SUB Square (a, b, c, d) PRINT PRINT " v4 v3" PRINT " "; d; " -----------"; c PRINT " | |" PRINT " | |" PRINT " | |" PRINT " | |" PRINT " | |" PRINT " "; a; " -----------"; b PRINT " v1 v2 " PRINT END SUB