'Program LEPLOT.BAS makes plot of largest LE for Smith, Thiffeault, Horton model vs. Vsw
'(c) 2000 by J. C. Sprott <http://sprott.physics.wisc.edu/>
'Compile with PowerBASIC for DOS 3.5 <http://www.powerbasic.com/>
'$error all on
SCREEN 12
defext a - z
shared a1, a2, b1, b2, b3, d1, f1, g1, g2, g3, vsw
randomize timer
line(0, 0) - (639, 479), ,b
for i%=1 to 9
    line(64*i%,479)-(64*i%,469)
    line(64*i%,1)-(64*i%,11)
    line(1,48*i%)-(11,48*i%)
    line(629,48*i%)-(639,48*i%)
next i%
pset(0,240)
call insolve(1, 1, 638, 478)
while inkey$ <> chr$(27): wend
end

SUB calculate(xp%, yp%, le!)
    h = 0.1     'Iteration step size
    dr = 1e-8   'Perturbation for LE calculation
    x = .06: y = .11: z = 43: w = .75
    xe = x + dr: ye = y: ze = z: we = w
    a1 =.247
    a2 =.391
    b1=10.8
    b2=0.0752
    b3=1.06
    d1=2200
    f1=2.47
    g1=1080
    g2=4
    g3=3.79
    vsw=.01+9.99*xp%/639
    iter& = 0: lsum = 0: ub% = 0
    for j& = 1 to 200000
    	call rk4(x, y, z, w, h)
        call rk4(xe, ye, ze, we, h)
        dx = xe - x: dy = ye - y: dz = ze - z: dw = we - w
        dsq = dx * dx + dy * dy + dz * dz + dw * dw
        rs = sqr(dsq) / dr
        if rs > 0 then
            xe = x + dx / rs: ye = y + dy / rs: ze = z + dz / rs: we = w + dw/rs
            if j& > 2000 then lsum = lsum + log(rs): incr iter&
        end if
        if abs(x) > 10000 then lsum=0: exit for
    next j&
    le! = lsum / (iter& * h)
    'locate 1,1: print le!
END SUB

def fnx(x,y,z,w) = a1*(vsw-y)+a2*(y-fnvi(y))
def fny(x,y,z,w) = b1*(x-fni1(y))-b2*sqr(z)-b3*y
def fnz(x,y,z,w) = y*y-sqr(w)*z*.5*(1+tanh(d1*(x-1)))
def fnw(x,y,z,w) = sqr(z)*y-w

def fnvi(y) = y+a2*(vsw-y)/f1
def fni1(y) = g3*g3*fnvi(y)^3/(2*g1*g1)+g2*fnvi(y)/g1+(g3*fnvi(y)^2/(2*g1*g1))*sqr(4*g1*g2+g3*g3*fnvi(y)^2)

SUB rk4 (x, y, z, w, h)
    'Advances (x, y, z, w) with a fourth order Runge-Kutta method
    'User must provide the functions fnx, fny, fnz, fnw, and a step size h
    k1x = h * fnx(x, y, z, w)
    k1y = h * fny(x, y, z, w)
    k1z = h * fnz(x, y, z, w)
    k1w = h * fnw(x, y, z, w)
    k2x = h * fnx(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z, w + .5 * k1w)
    k2y = h * fny(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z, w + .5 * k1w)
    k2z = h * fnz(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z, w + .5 * k1w)
    k2w = h * fnw(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z, w + .5 * k1w)
    k3x = h * fnx(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z, w + .5 * k2w)
    k3y = h * fny(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z, w + .5 * k2w)
    k3z = h * fnz(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z, w + .5 * k2w)
    k3w = h * fnw(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z, w + .5 * k2w)
    k4x = h * fnx(x + k3x, y + k3y, z + k3z, w + k3w)
    k4y = h * fny(x + k3x, y + k3y, z + k3z, w + k3w)
    k4z = h * fnz(x + k3x, y + k3y, z + k3z, w + k3w)
    k4w = h * fnw(x + k3x, y + k3y, z + k3z, w + k3w)
    x = x + (k1x + 2 * (k2x + k3x) + k4x) / 6
    y = y + (k1y + 2 * (k2y + k3y) + k4y) / 6
    z = z + (k1z + 2 * (k2z + k3z) + k4z) / 6
    w = w + (k1w + 2 * (k2w + k3w) + k4w) / 6
END SUB

SUB insolve (x1%, y1%, x2%, y2%)
    'Makes a plot of LE vs Vsw
    for x%=x1% to x2% 'step 10
    	if inkey$=chr$(27) then end
    	call calculate (x%,y%,le!)
        line-(x%,240-480*le!)
    next x%
END SUB

FUNCTION tanh (x)
    'Returns hyperbolic tangent of x
    tanh = 1 - 2 / (EXP(2 * x) + 1)
END FUNCTION