REM
REM SSAFR18.bas Large loop model including (dlim_0, dlim_c, dlim_p).
REM Peter Lohmander 170708_2200
REM Software references:
REM tutorial: https://www.youtube.com/watch?v=GjrSxMUb9F8
REM http://www.qb64.net/wiki/index.php/COLOR
DIM x(1000), y(1000), d(1000), a(1000), c(1000), dist(1000, 1000)
DIM h(1000), vol(1000), p(1000), net(1000), csmall(1000), qual(1000)
DIM pdev(100, 200), height(1000), dadt(1000)
CLS
OPEN "a_OutSSAFR18.txt" FOR OUTPUT AS #2
PRINT #2, "Results from SSAFR18.bas by Peter Lohmander 170708"
PRINT "Results from SSAFR18.bas by Peter Lohmander 170708"
PRINT #2, " r SigmaK dlim_0 dlim_c dlim_q dlim_p avHARV TotalPV"
PRINT " r SigmaK dlim_0 dlim_c dlim_q dlim_p avHARV TotalPV"
REM parameters
PI = 3.1415926536
dmax = 10
hacorr = 10000 / (PI * (dmax / 2) ^ 2)
formnumb = 0.5
setupc = 1000
r = 0.03
TMAX = 200
REM Harvest parameter definitions
harvt = 5
REM Generation of random market price deviations
randseed = 1
RANDOMIZE randseed
Numser = 10
PEQU = 503.3
SigmaP = 59.03
FOR nums = 1 TO Numser
Price_t = PEQU
pdev(nums, 0) = 0
FOR t = 1 TO 200
REM Generation of Random Number N(0,1):
epsilon = 0
FOR ii = 1 TO 12
epsilon = epsilon + RND
NEXT ii
epsilon = (epsilon - 6)
Price_t = 235.9 + 0.5313 * Price_t + SigmaP * epsilon
pdev(nums, t) = Price_t - PEQU
NEXT t
NEXT nums
REM Here the loop starts with alternative values of
REM some important exogenous parameters.
REM FOR r = 0.02 TO 0.041 STEP 0.01
r = 0.03
REM FOR SigmaK = 0.5 TO 1.51 STEP 0.5
SigmaK = 1
REM Here is the loop where all total system simulations are made
REM with different control function parameter values
FOR dlim_0 = 0.20 TO 0.51 STEP .050
FOR dlim_c = 0.0000 TO -0.0036 STEP -0.0005
REM FOR dlim_q = 0.00 TO 0.061 STEP 0.02
dlim_q = 0.0600
FOR dlim_p = -0.0070 TO 0.0001 STEP 0.001
Detail = 0
sum_Totharv = 0
sum_Totpresv = 0
FOR nums = 1 TO Numser
GOSUB Subroutine_SDOC
sum_Totharv = sum_Totharv + Totharv
sum_Totpresv = sum_Totpresv + Totpresv
NEXT nums
PRINT #2, USING "###.####"; r;
PRINT #2, USING "###.####"; SigmaK;
PRINT #2, USING "###.####"; dlim_0;
PRINT #2, USING "###.####"; dlim_c;
PRINT #2, USING "###.####"; dlim_q;
PRINT #2, USING "###.####"; dlim_p;
PRINT #2, USING "####.####"; sum_Totharv / 200 / Numser;
PRINT #2, USING "#######.###"; sum_Totpresv / Numser
PRINT USING "###.####"; r;
PRINT USING "###.####"; SigmaK;
PRINT USING "###.####"; dlim_0;
PRINT USING "###.####"; dlim_c;
PRINT USING "###.####"; dlim_q;
PRINT USING "###.####"; dlim_p;
PRINT USING "####.####"; sum_Totharv / 200 / Numser;
PRINT USING "#######.###"; sum_Totpresv / Numser
NEXT dlim_p
REM NEXT dlim_q
NEXT dlim_c
NEXT dlim_0
REM NEXT SigmaK
REM NEXT r
CLOSE #2
END
REM Here the main program ends.
REM Here, the subroutine starts, where all of the full system simulations are made
Subroutine_SDOC:
IF Detail = 1 THEN SCREEN _NEWIMAGE(800, 800, 256)
IF Detail = 1 THEN COLOR 30, 15
REM Initial conditions via diameter frequences
REM The following data are relevant to the Romperod forest area.
REM In Rompedod, there are 825 trees per hectare. Here, however, we
REM assume that there are 1000 trees per hectare, but with the same
REM relative frequency distribution.
randseed = 1
RANDOMIZE randseed
FOR N = 1 TO 1000
x(N) = 100 * RND
y(N) = 100 * RND
dialim = 0
random1 = RND
IF random1 > 0.3266 THEN dialim = 0.1
IF random1 > 0.5859 THEN dialim = 0.2
IF random1 > 0.7508 THEN dialim = 0.3
IF random1 > 0.8788 THEN dialim = 0.4
IF random1 > 0.9663 THEN dialim = 0.5
IF random1 > 0.9966 THEN dialim = 0.6
d(N) = dialim + RND * 0.1
a(N) = PI * (d(N) / 2) ^ 2
qual(N) = 0
IF RND > .3 THEN qual(N) = 1
NEXT N
FOR n1 = 1 TO 1000
FOR n2 = 1 TO 1000
dist(n1, n2) = ((x(n1) - x(n2)) ^ 2 + (y(n1) - y(n2)) ^ 2) ^ .5
NEXT n2
NEXT n1
Totharv = 0
Totpresv = 0
FOR t = 1 TO TMAX
IF Detail = 1 THEN CLS
IF Detail = 1 THEN COLOR 3: PRINT ""
IF Detail = 1 THEN COLOR 3: PRINT " t = "; t
harvper = 0
IF (INT(t / harvt) = t / harvt) THEN harvper = 1
IF harvper = 1 AND Detail = 1 THEN PRINT " This year, harvest is possible."
REM Calculation of competition for individual n1.
REM Only trees closer than dmax are considered.
REM
FOR n1 = 1 TO 1000
sumAT = 0
sumS = 0
FOR n2 = 1 TO 1000
IF n2 = n1 THEN GOTO 100
IF dist(n1, n2) > dmax THEN GOTO 100
sumAT = sumAT + (a(n2) / a(n1)) ^ 0.1515 * a(n2) * EXP(-dist(n1, n2) / 5.630) ^ 2
sumS = sumS + (a(n2) / a(n1)) ^ 0.1515 * a(n2) * EXP(-dist(n1, n2) / 5.630) ^ 2
100 NEXT n2
c(n1) = 83.16 * sumAT ^ 1.1768 + 55.96 * sumS ^ 1.1768
NEXT n1
REM Harvest decisions
IF harvper = 0 THEN GOTO 200
FOR n1 = 1 TO 1000
h(n1) = 0
IF d(n1) > (dlim_0 + dlim_c * c(n1) + dlim_q * qual(n1) + dlim_p * SigmaK * pdev(nums, t)) THEN h(n1) = 1
NEXT n1
REM Harvest volumes and values etc
Hvolume = 0
Hnet = 0
Hpresv = 0
FOR n1 = 1 TO 1000
IF h(n1) = 0 THEN GOTO 500
diam_mm = d(n1) * 1000
IF diam_mm > 674 THEN diam_mm = 674
hojd_dm = 13 + 1.4743 * diam_mm - 0.037864 * diam_mm ^ 1.5
height(n1) = hojd_dm / 10
vol(n1) = (PI * (d(n1) / 2) ^ 2) * height(n1) * formnumb
Hvolume = Hvolume + vol(n1)
p(n1) = (PEQU + SigmaK * pdev(nums, t) + 150 * (qual(n1) - 0.5) + 2 * (d(n1) - 40)) * (0.8 * 1.0 + 0.2 * 0.6) - 206
net(n1) = p(n1) * vol(n1)
Hnet = Hnet + net(n1)
500 REM
NEXT n1
Totharv = Totharv + Hvolume
Hpresv = EXP(-r * t) * (-setupc + Hnet)
Totpresv = Totpresv + Hpresv
IF Detail = 1 THEN PRINT " Harvest volume this year = "; Hvolume;
IF Detail = 1 THEN PRINT " Present value of activitites this year = "; Hpresv
IF Detail = 1 THEN PRINT " Total harvest volume= "; Totharv;
IF Detail = 1 THEN PRINT " Total present value = "; Totpresv
REM Consequences after harvesting:
REM In case n1 has been harvested, the individual n1
REM is removed and a individual with random coordinates
REM appears in the list with small size
FOR n1 = 1 TO 1000
IF h(n1) = 0 THEN GOTO 300
x(n1) = 100 * RND
y(n1) = 100 * RND
d(n1) = 0.1 * RND
a(n1) = PI * (d(n1) / 2) ^ 2
qual(n1) = 0
IF RND > .5 THEN qual(n1) = 1
300 REM
NEXT n1
REM Calculation of new distances between n1 and n2
REM in case at least one of them has been harvested.
FOR n1 = 1 TO 1000
FOR n2 = 1 TO 1000
newdist = h(n1) + h(n2)
IF newdist < 1 THEN GOTO 400
dist(n1, n2) = ((x(n1) - x(n2)) ^ 2 + (y(n1) - y(n2)) ^ 2) ^ .5
400 REM
NEXT n2
NEXT n1
200 REM
REM Calculation of growth of all individuals n.
REM first, the basal area increment dadt(n) (m2/year), is calculated.
REM Then, the basal area after growth, a(n) (m2) is calculated.
REM Finally, the diameter after growth, d(n) (m), is calculated.
FOR N = 1 TO 1000
dadt(N) = (88.50 + 73.03) / 2 * a(N) ^ 0.5 + (-154.71 - 85.33) / 2 * a(N) ^ 1.5
dadt(N) = dadt(N) - c(N) * a(N) ^ 0.5
dadt(N) = dadt(N) / 10000
IF dadt(N) < 0 THEN dadt(N) = 0
a(N) = a(N) + dadt(N)
d(N) = (4 * a(N) / PI) ^ 0.5
NEXT N
IF Detail = 0 THEN GOTO 600
FOR N = 1 TO 1000
CIRCLE (INT(6 * x(N) + 100), INT(6 * y(N) + 100)), INT(100 * d(N) / 5), 0
NEXT N
LINE (100, 100)-(700, 100), 3
LINE (100, 100)-(100, 700), 3
LINE (700, 100)-(700, 700), 3
LINE (100, 700)-(700, 700), 3
SLEEP
600 REM
NEXT t
RETURN