'Program WEALTH.BAS (c) 2012 by J. C. Sprott
'Calculates the distribution of wealth in a society subject to certain assumptions
'Use PowerBasic Console Compiler
#DEBUG ERROR ON
DEFEXT a-z

FUNCTION PBMAIN()
w&=500      'width of graph
h&=500      'height of graph
p&=1e6      'population size
z&=1e5      'initial wealth
y&=20       'number of years
r&=20       'range of rates (percent)
m&=1e6      'maximum wealth
d&=2        'decades to plot
CONSOLE SET LOC w&+6,0
CLS
GRAPHIC WINDOW EXE.NAME$, 0, 0, w&, h& TO hWin???
GRAPHIC ATTACH hWin???, 0 ', REDRAW
GRAPHIC COLOR %BLACK, %WHITE
GRAPHIC CLEAR
GRAPHIC BOX(0,0)-(w&,h&)
CONSOLE SET FOCUS
RANDOMIZE TIMER
DIM wealth(p&),rate(p&),wbin(w&)
FOR i&=1 TO p&
    wealth(i&)=z&
NEXT i&
FOR i&=1 TO p&
    'rate(i&)=0.01##*r&*(RND-0.5##)  'uniform random
    rate(i&)=0.005##*r&*gauss       'Gaussian random
    FOR year&=1 TO y&
        'rate=0.01##*r&*(rnd-0.5##)      'different rate every year
        rate=rate(i&)                   'same rate every year
        'rate=0.5##*rate(i&)+0.5##*0.01##*r&*(RND-0.5##) 'mixed rate
        'rate=0.5##*rate(i&)+0.5##*0.005##*r&*gauss      'mixed rate Gaussian
        wealth(i&)=(1##+rate)*wealth(i&)
    NEXT year&
    wb&=INT(w&+w&*LOG10(wealth(i&)/m&)/d&)
    IF wb&>0 AND wb&<=w& THEN
        wbin(wb&)=wbin(wb&)+10^(-d&*wb&/w&)
        wbm&=MAX(wbm&,wbin(wb&))
    END IF
NEXT i&
FOR x&=1 TO w&
    IF wbin(x&)=0 THEN
        y&=h&
    ELSE
        y&=-h&*LOG10(wbin(x&)/wbm&)/d&
    END IF
    IF x&=1 THEN GRAPHIC SET PIXEL(x&,y&) ELSE GRAPHIC LINE-(x&,y&)
    PRINT x&,y&
NEXT w&
GRAPHIC SAVE"wealth4.bmp"
WAITKEY$
END FUNCTION

FUNCTION gauss
    'Returns a normally (Gaussian) distributed random deviate
    'with mean of zero and standard deviation of 1.0
    DO
        V1 = 2## * RND - 1##
        V2 = 2## * RND - 1##
        R = V1 * V1 + V2 * V2
    LOOP WHILE R >= 1## OR R = 0##
    GAUSS = V2 * SQR(-2## * LOG(R) / R)
END FUNCTION