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