DECLARE SUB pager ()
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                                                               
1200  REM                                                               
1300  REM  description:                                                 
1400  REM                                                               
1500  REM       this program implements the strategy developed for the  
1600  REM  nuclear decay problem and applies it to ecology.  The example
1700  REM  is two species with different growth rates and interactions. 
1800  REM                                                               
1900  REM                                                               
2000  REM  instructions:                                                
2100  REM                                                               
2200  REM       this program is taken from "using computers in          
2300  REM  physics" by john r. merrill and should be used in conjunction
2400  REM  with that text.  The programs are distributed by project     
2500  REM  conduit.                                                   
REM extensively modified by
REM GE Copeland
REM Dept of Physics
REM Old Dominion University
REM Norfolk, VA 23504
REM 3/28/2001
REM


2505  PRINT CHR$(12);
2600 PRINT , "Title:ecology"
2700 PRINT , "     An ecology simulation program"
2800 PRINT
2900 PRINT
3000  REM                                                               
3100 PRINT "The program 'Ecology' implements the strategy as developed"
3200 PRINT "for the nuclear-decay problem (see programs U234.bas and/"
3300 PRINT "or Nucdec.bas for an introduction to general decay"
3400 PRINT "concepts), but applied to an ecological situation. Since"
3500 PRINT "ecological equations for populations of species living in regions which do not communicate"
3600 PRINT "with the rest of the world are quite similar to the rate"
3700 PRINT "equations that govern nuclear decay problems."
3800 PRINT
3900 CALL pager
4000 PRINT "The basic mathematical equations for two interacting"
4100 PRINT "species are the Lotka-Volterra equations:"
4200 PRINT
4300 PRINT "      (dn1/dt)=r1*n1*(k1-n1-j1*n2)/k1"
4400 PRINT "      (dn2/dt)=r2*n2*(k2-n2-j2*n1)/k2"
4500 PRINT
4600 PRINT "the ri's are the 'biotic potentials' (unrestricted"
4700 PRINT "growth-rate constants). The ki's are the 'carrying"
4800 PRINT "capacities' (max. Number of each species the region can"
4900 PRINT "support). The j's are the coefficients of species"
5000 PRINT "interaction, which are positive for inhibitory influences"
5100 PRINT "(i.e. one species eats the other) and negative for effects"
5200 PRINT "which increase numbers."
5300 CALL pager
5400 PRINT
5500 PRINT "Rather than do a complete derivation, let's consider"
5600 PRINT "what reasonable terms in an ecological equation might"
5700 PRINT "look like. First, consider a region with a large "
5800 PRINT "carrying capacity and no species competition. The rate"
5900 PRINT "at which the species will increase is just proportional"
6000 PRINT "to the number of mature adults. For simplicity, assume"
6100 PRINT "maturation is instantaneous (number of adults=total"
6200 PRINT "number of individuals). Then:"
6300 PRINT
6400 PRINT "        (dn1/dt)=c1*n1, and  c1 >0"
6500 PRINT
6600 PRINT "note: if we rearranged this, dn1/n1=cdt, and integrated"
6700 PRINT "      (Limits of 1,n,0,t) we would obtain our basic nuclear"
6800 PRINT "      decay equation, excepting c was negative in that"
6900 PRINT "      situation."
7000 CALL pager
7200 PRINT "As a next approximation let's assume a finite carrying"
7300 PRINT "capacity, whose max number is m1. You must reflect this"
7400 PRINT "in your equation. Hence"
7500 PRINT
7600 PRINT "      (dn1/dt)=c1*n1*(1-n1/m1)"
7700 PRINT
7800 PRINT "for n1<<m1 the growth is unaffected. However, as n1"
7900 PRINT "approaches m1 the rate goes to zero. This gives us the"
8000 PRINT "end values we need, but in between we don't approach"
8100 PRINT "nature realistically (too simplified). Thus you should"
8110 PRINT "consider the curve-rate equation as the beginning of a"
8120 PRINT "series expansion."
8200 PRINT
8300 PRINT "A third approximation should include species competition."
8400 PRINT "for simplicity, one assumes the interactions are propor-"
8500 PRINT "tional to the numbers of individuals in each species"
8600 PRINT "(unlike decay problems). These interactions should depend"
8700 PRINT "on the species meeting (a function of how many of each"
8800 CALL pager
8900 PRINT "species there are). If the interaction term is written"
9000 PRINT "as a1n1n2, then the system of equations for two com-"
9100 PRINT "petitive species is:"
9200 PRINT
9300 PRINT "       (dn1/dt)=c1*n1*(1-n1/m1)+a1*n1*n2"
9400 PRINT "       (dn2/dt)=c2*n2*(1-n2/m2)+a2*n2*n1"
9500 PRINT
9600 PRINT "where ci and mi are positive, but ai can have either"
9700 PRINT "sign. This form of the equation is merely that of"
9800 PRINT "Lotka Volterra if you write ci=ri, mi=ki and"
9900 PRINT "ai=-ji ri/ki. In the ecology case (unlike the nuclear"
10000 PRINT "decay case) there is no necessary relationship between"
10100 PRINT "a1 and a2 (i.e. it may not help the lion as much to eat"
10200 PRINT "a man as it hurts the man to be eaten)."
10400 CALL pager
10500 PRINT "By now we are sure you're impatient to be into the program,"
10600 PRINT "just a few more comments and we will begin the calculation. "
10800 PRINT
10900 PRINT "This program will assume you start with two members of"
11000 PRINT "each species with growth rates of three and one per year."
11100 PRINT "The max carrying capacities are 1e+5 and 1e+2,respectively."
11200 PRINT "But the interaction of species 2 on species 1 (a1) will be"
11300 PRINT "left for you to select. For simplicity, the interaction"
11400 PRINT "of species 1 on 2 is arbitrarily set at zero. Also, the"
11500 PRINT "program assumes that once a species is depleted (n<1),"
11600 PRINT "it does not recover."
12000  REM  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *   
12100 DIM t9(300), c1(300), c2(300)
12200 n9 = 0
12300  READ n1, n2
12400  DATA 2,2                                                          
12500 READ r1, m1, r2, m2, a2
12600 DATA 5,100000,1,100,0                                              
12700 PRINT "Do you want to pick a value for a1 (yes or no)";
12800 INPUT b$
12900 IF b$ = "no" THEN GOTO 14100
13000 IF b$ = "yes" THEN GOTO 13400
13100 REM shouldn't be here                                              
13200 PRINT "I don't understand."
13300 GOTO 12700
13400 PRINT "You must now choose a value between -.3 and +.3 (the limits"
13500 PRINT "are designed for convenance and simplicity)."
13600 INPUT f1
13700 IF f1 > .3 THEN GOTO 13400
13800 IF f1 < -.3 THEN GOTO 13400
13900  a1 = f1
14000 GOTO 14300
14100 PRINT "We set a1=-.1."
14200  a1 = -.1
14300   t = 0
14400 PRINT "Do you want to see numerical values (yes or no)";
14500 INPUT z$
14600 IF z$ = "yes" THEN GOTO 15100
14700 IF z$ = "no" THEN GOTO 15100
14800 REM shouldn't be here                                              
14900 PRINT "i don't understand"
15000 GOTO 14400
15100   d = .01
ik = 0
15200   t1 = .499
15250 IF z$ = "no" THEN GOTO 15600
15300 PRINT
15400 PRINT "Time/yr", "Species i", "100*Species ii"
15500 PRINT
15600   l1 = n1 + (r1 * n1 * (1 - n1 / m1) + a1 * n1 * n2) * d
15700   l2 = n2 + (r2 * n2 * (1 - n2 / m2) + a2 * n2 * n1) * d
15800   l1 = (l1 + n1) / 2
15900   l2 = (l2 + n2) / 2
16000   n1 = n1 + (r1 * l1 * (1 - l1 / m1) + a1 * l1 * l2) * d
16100  IF n1 > 1 THEN GOTO 16300
16200   n1 = 0
16300   n2 = n2 + (r2 * l2 * (1 - l2 / m2) + a2 * l2 * l1) * d
16400  IF n2 > 1 THEN GOTO 16600
16500   n2 = 0
16600   t = t + d
16700  IF t <= t1 THEN GOTO 17600
16800 IF z$ = "yes" THEN GOTO 17000
16900 IF z$ = "no" THEN GOTO 17100
17000
      IF ik = 10 THEN
      PRINT t, n1, 100 * n2
      ik = 0
      ELSE
      ik = ik + 1
      END IF

17100 n9 = n9 + 1
17200 t9(n9) = t
17300 c1(n9) = n1
17400 c2(n9) = n2
17500   t1 = t1 + .05
17600  IF t <= 12.5 THEN GOTO 15600
17700 PRINT "Would you like a graph of the data (yes or no)";
17800 INPUT h$
17900 IF LEFT$(h$, 1) = "y" OR LEFT$(h$, 1) = "Y" THEN GOTO 18900
18000 IF h$ = "no" THEN STOP
18100 REM shouldn't be here                                              
18200 PRINT "i don't understand"
18300 GOTO 17700
18900 REM output to plttek.dat file for graphics plotting                
19000 OPEN "plttek.dat" FOR OUTPUT AS #2
19200 PRINT #2, 1
REM !# of sheets                                        
19300 PRINT #2, 2
REM !# of graphs on first sheet                         
19400 PRINT #2, n9
REM !# of data pairs per graph                          
19500 PRINT #2, "      Biological units"
19600 PRINT #2, "     time(years)--->"
19700 PRINT #2, "     Number of biological units as a"
REM !title                      
19800  PRINT #2, "      function of time"
19900 FOR I = 1 TO n9
20000 PRINT #2, t9(I), c1(I)
20100 NEXT I
20200 FOR I = 1 TO n9
20300 PRINT #2, t9(I), 100 * c2(I)
20400 NEXT I
CLOSE #2

20500 PRINT
20600 PRINT "Now run WPLOT."
CHAIN "Tek2wplt.exe"
20700 STOP
22600 STOP
22700 PRINT "this completes our program."
22800  END

SUB pager
PRINT , "Push return to continue.";
INPUT r$
PRINT CHR$(12);
END SUB