REM REM SpaFor.bas REM Peter Lohmander 160707 REM tutorial: https://www.youtube.com/watch?v=GjrSxMUb9F8 REM http://www.qb64.net/wiki/index.php/COLOR REM Some equations from Schutz article REM i(k)=b0+b1*LN(dk)+b2(C(k)) ^3 REM i(k)= diameter increment per year REM diameter before increment REM total basal area per hectare of trees with diameter larger than dk REM b0 = 1.506969 REM b1 = 0.94255 REM b2 = -0.000183455 DIM x(1000), y(1000), d(1000), a(1000), i(1000), c(1000), dist(1000, 1000) DIM h(1000), vol(1000), p(1000), net(1000), csmall(1000), qual(1000) DIM pdev(100, 200) CLS OPEN "a_Out26.txt" FOR OUTPUT AS #2 PRINT #2, " dlim_0 dlim_c dlim_q dlim_p avHARV TotalPV" REM parameters PI = 3.1415926536 b0 = 1.506969 b1 = 0.94255 b2 = -0.000183455 dmax = 10 hacorr = 10000 / (PI * (dmax / 2) ^ 2) hdratio = 50 formnumb = 0.5 setupc = 100 p0 = -20 p1 = 100 p2 = 40 pqual = 15 r = 0.03 TMAX = 200 REM Harvest parameter definitions dlim = 0.40 harvt = 5 REM Generation of random market price deviations randseed = 1 RANDOMIZE randseed Numser = 10 rmp = 10 FOR t = 0 TO 200 FOR nums = 1 TO Numser pdev(nums, t) = rmp * (2 * RND - 1) NEXT nums NEXT t REM Here is the loop where all total system simulations are made REM with different control function parameter values FOR dlim_0 = 0.40 TO 0.81 STEP .025 dlim_c = -0.003 REM FOR dlim_c = -.007 TO 0 STEP 0.002 dlim_q = 0.020 REM FOR dlim_q = 0.02 TO 0.04 STEP 0.01 FOR dlim_p = -0.040 TO 0.001 STEP 0.005 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 " dlim_0= "; dlim_0; " dlim_c= "; dlim_c; PRINT " dlim_q= "; dlim_q; " dlim_p= "; dlim_p; PRINT " PV= "; sum_Totpresv / Numser 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 NEXT dlim_p REM NEXT dlim_q REM NEXT dlim_c NEXT dlim_0 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 FOR N = 1 TO 1000 x(N) = 100 * RND y(N) = 100 * RND d(N) = 0.43 * RND a(N) = PI * (d(N) / 2) ^ 2 qual(N) = 0 IF RND > .5 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 = -2 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 according to the Schutz function but where REM only trees closer than dmax are considered. REM FOR n1 = 1 TO 1000 c(n1) = 0 FOR n2 = 1 TO 1000 IF n2 = n1 THEN GOTO 100 IF dist(n1, n2) > dmax THEN GOTO 100 IF d(n2) < d(n1) THEN GOTO 100 c(n1) = c(n1) + a(n2) 100 NEXT n2 c(n1) = c(n1) * hacorr NEXT n1 REM Calculation of competition for individual n1 REM if only smaller trees than the individual tree and REM only trees closer than dmax are considered. REM FOR n1 = 1 TO 1000 csmall(n1) = 0 FOR n2 = 1 TO 1000 IF n2 = n1 THEN GOTO 101 IF dist(n1, n2) > dmax THEN GOTO 101 IF d(n2) > d(n1) THEN GOTO 101 csmall(n1) = csmall(n1) + a(n2) 101 NEXT n2 csmall(n1) = csmall(n1) * hacorr 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) + csmall(n1)) + dlim_q * qual(n1) + dlim_p * 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 h(n1) = hdratio * d(n1) vol(n1) = ((d(n1) / 2) ^ 2) * h(n1) * formnumb Hvolume = Hvolume + vol(n1) p(n1) = p0 + p1 * d(n1) + pqual * qual(n1) IF p(n1) > p2 THEN p(n1) = p2 p(n1) = p(n1) + pdev(nums, t) 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 according REM to the competition distance adjusted Schutz function REM i(k)=b0+b1*LN(dk)+b2(C(k)) ^3 REM Important modification: REM Note that (csmall(.)/2) has been included by PL in the function! FOR N = 1 TO 1000 i(N) = b0 + b1 * LOG(1000 * d(N)) + b2 * (c(N) + csmall(N) / 2) ^ 3 IF i(N) < 0 THEN i(N) = 0 d(N) = d(N) + i(N) / 1000 a(N) = PI * (d(N) / 2) ^ 2 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