!   CASSINI.TB

!   This TrueBasic program introduces the basic framework of classical 
!   astronomy, then shows (schematically) a cross section of the basilica 
!   of San Petronio and a set of rays from the sun passing through the 
!   hole inthe roof of the basilica and tracing an arc on the floor of
!   the vault.

!   INITIAL DECLARATIONS

    dim dn(3)                           !  Direction numbers
    dim io(3), jo(3), ko(3)             !  Old basis vectors
    dim in(3), jn(3), kn(3)             !  New basis vectors
    dim xsun(3)                         !  Start position for sun
    dim ss(2), s1(2), s2(2), s3(2)      !  Variable vectors
    dim uu(3), vv(3), ww(3)             !  Variable vectors
    dim w1(3), w2(3), w3(3)             !  Variable vectors
    dim xx(3), yy(3), zz(3)             !  Variable vectors
    dim xs(3), ys(3), zs(3)             !  Variable vectors
    dim x1(3), x2(3), x3(3)             !  Variable vectors
    dim y1(3), y2(3), y3(3)             !  Variable vectors
    dim z1(3), z2(3), z3(3)             !  Variable vectors
    dim zp(3), cp(3), ep(3)             !  Variable vectors
    dim aa(3,3), bb(3,3), cc(3,3)       !  Variable matrices
    dim dd(3,3), ee(3,3)                !  Variable matrices
    dim hh(3,3), kk(3,3), mm(3,3)       !  Fixed matrices
    dim pp(3,3), qq(3,3), rr(3,3)       !  Variable matrices
    dim id(3,3)                         !  Identity matrix
    dim eq(3,3), ec(3,3), nn(3,3)       !  Transformation matrices

!   INPUT OPTIONS

    option angle degrees
    let io(1) = 1
    let jo(2) = 1
    let ko(3) = 1
    mat zp = ko
    let rad = 0.01
    for jj = 1 to 3
        let id(jj,jj) = 1
    next jj
!   input prompt "Set direction numbers for eye-point ": dn(1), dn(2), dn(3)
!   input prompt "Set distance for eye-point ": eye
    let dn(1) = 3
    let dn(2) = -3
    let dn(3) = 1
    let eye = 5
!   let eyeinv = 1/eye
    let lt = 44.5   !  The latitude of Bologna.
    let epsilon = 23.5    !   The obliquity of the ecliptic.

!   PRELIMINARY COMPUTATIONS

    call unit(dn(),kn())
    call cross(ko(),kn(),zz())
    call unit(zz(),in())
    call cross(kn(),in(),jn())
    for jj = 1 to 3
        let nn(1,jj) = in(jj)
        let nn(2,jj) = jn(jj)
        let nn(3,jj) = kn(jj)
    next jj
    !  We have yy = nn*xx, where xx are the coordinates of a vector relative
    !  to the standard basis io,jo,ko and yy are the coordinates of that vector 
    !  relative to the basis of perspective in,jn,kn.
    
!   WINDOW SETTINGS

    ask pixels px, py
    let wx = (px/py)*1.2
    let wy = 1.2
    let vx = (px/py)*1.1
    let vy = 1.1
    set window -wx, wx, -wy, wy
    set color 255    !  black
    box area -wx, wx, -wy, wy
    set color 0    !  white
    box lines -vx, vx, -vy, vy
    !  box circle -1,1,-1,1

!   BOXES

    set color 255   !  black
    box circle -3*rad,3*rad,-3*rad,3*rad
    flood 0,0
    box keep -3*rad,3*rad,-3*rad,3*rad in blk$    !  black box
    set color 155   !  green
    box circle -3*rad,3*rad,-3*rad,3*rad
    flood 0,0
    box keep -3*rad,3*rad,-3*rad,3*rad in earth$    !  earth box
    box show blk$ at -3*rad,-3*rad  
    box circle -rad,rad,-rad,rad
    flood 0,0
    box keep -rad,rad,-rad,rad in zpole$    !  zenith box
    box show blk$ at -3*rad,-3*rad
    set color 17   !  gold
    box circle -2*rad,2*rad,-2*rad,2*rad
    flood 0,0
    box keep -2*rad,2*rad,-2*rad,2*rad in sun1$    !  sun box
    box show blk$ at -3*rad,-3*rad
    box circle -rad,rad,-rad,rad
    flood 0,0
    box keep -rad,rad,-rad,rad in sun2$    !  sun spot
    box show blk$ at -3*rad,-3*rad
    set color 198   !  blue
    box circle -rad,rad,-rad,rad
    flood 0,0
    box keep -rad,rad,-rad,rad in cpole$    !  celestial pole box
    box show blk$ at -3*rad,-3*rad
    set color 35   !  red
    box circle -rad,rad,-rad,rad
    flood 0,0
    box keep -rad,rad,-rad,rad in epole$    !  ecliptic pole box
    box show blk$ at -3*rad,-3*rad

!   FRAMEWORK

    call polar(lt,0,zz())
    call cross(jo(),zz(),xx())
    for jj = 1 to 3
        let eq(jj,1) = xx(jj) 
        let eq(jj,2) = jo(jj) 
        let eq(jj,3) = zz(jj) 
    next jj    
    set color mix(249) 0.2,0.2,0.2
    for th = -75 to 75 step 15
    if th = 0 then
        set color 198    !  blue
    else
        set color 249    !   dark gray
    end if
    for ph = 0 to 361 step 3
        call polar(th,ph,xx())
        mat yy = eq*xx
        mat zz = nn*yy
        call persp(zz(),ss())
        plot lines: ss(1),ss(2);
    next ph
        plot points: ss(1), ss(2)
    next th

    set color 249    !  dark gray
    for ph = 0 to 330 step 45
    for th = 0 to 75 step 3
        call polar(th,ph,xx())
        mat yy=eq*xx
        mat zz = nn*yy
        call persp(zz(),ss())
        plot lines: ss(1),ss(2);
    next th
        plot points: ss(1), ss(2)
    next ph

    set color 155    !  green
    for ph = 0 to 361 step 3
        call polar(0,ph,yy())
        mat zz = nn*yy
        call persp(zz(),ss())
        plot lines: ss(1),ss(2);
    next ph
        plot points: ss(1), ss(2)

    set color 0    !  white
    for th = 0 to 90 step 3
        call polar(th,0,yy())
        mat zz = nn*yy
        call persp(zz(),ss())
        plot lines: ss(1),ss(2);
    next th
        plot points: ss(1), ss(2)

    set color 0    !  white
    for th = 0 to 90 step 3
        call polar(th,180,yy())
        mat zz = nn*yy
        call persp(zz(),ss())
        plot lines: ss(1),ss(2);
    next th
        plot points: ss(1), ss(2)

    set color 35  !  red
    call polar(lt+epsilon,0,zz())
    call cross(jo(),zz(),xx())
    for jj = 1 to 3
        let ec(jj,1) = xx(jj) 
        let ec(jj,2) = jo(jj) 
        let ec(jj,3) = zz(jj) 
    next jj
    for ph = 0 to 361 step 3
        call polar(0,ph,xx())
        mat yy = ec*xx
        mat zz = nn*yy
        call persp(zz(),ss())
        plot lines: ss(1),ss(2);
    next ph
        plot points: ss(1), ss(2)
   
    set color 155   !  green
    mat vv = nn*zp
    call persp(vv(),s1())
    box show zpole$ at s1(1)-rad,s1(2)-rad
        plot text, at s1(1)-0.01, s1(2)+0.035: "Z"
    plot lines: 0,0; s1(1),s1(2)

    set color 155   !  green
    call polar(0,0,yy())
    mat zz = nn*yy
    call persp(zz(),ss())
    box show zpole$ at ss(1)-rad,ss(2)-rad

    set color 198   !  blue
    call polar(lt,0,cp())
    mat zz = nn*cp
    call persp(zz(),ss())
    box show cpole$ at ss(1)-rad,ss(2)-rad
        plot text, at ss(1)+0.015, ss(2)+0.03: "P"
    plot lines: 0,0; ss(1),ss(2)

    set color 35   !  red
    call polar(lt+23.5,0,ep())
    mat zz = nn*ep
    call persp(zz(),ss())
    box show epole$ at ss(1)-rad,ss(2)-rad
        plot text, at ss(1)-0.00, ss(2)+0.04: "Q"
    plot lines: 0,0; ss(1),ss(2)

    set color 155    !  green
    let th = 0
    for ph = 0 to 270 step 90
        call polar(th,ph,yy())
        mat zz = nn*yy
        call persp(zz(),ss())
        plot lines: 0,0; ss(1),ss(2)
        if ph =   0 then plot text, at ss(1)+0.01, ss(2)-0.06: "N"
        if ph =  90 then plot text, at ss(1)+0.02, ss(2)+0.02: "W"
        if ph = 180 then plot text, at ss(1)+0.025, ss(2)+0.03: "S"
        if ph = 270 then plot text, at ss(1)-0.035, ss(2)-0.06: "E"
    next ph

    set color 17   !  gold
    call polar(0,270,xx())        
    mat xsun = eq*xx        !  The sun begins at the spring equinox.
    mat zs = nn*xsun
    call persp(zs(),ss())

    box show earth$ at -3*rad,-3*rad

!   SAN PETRONIO (CROSS SECTION)

    set color 155    !   green
    let tc = 0.1
    let ta = 0.35
    let x1(1) = -ta
    let x1(2) = 0
    let x1(3) = 0
    let x2(1) = -ta
    let x2(2) = 0
    let x2(3) = tc
    mat y1 = nn*x1
    mat y2 = nn*x2
    call persp(y1(),s1())
    call persp(y2(),s2())
    plot lines: s1(1),s1(2); s2(1),s2(2)
    let x1(1) = -ta
    let x1(2) = 0
    let x1(3) = tc
    let x2(1) = ta
    let x2(2) = 0
    let x2(3) = tc
    mat y1 = nn*x1
    mat y2 = nn*x2
    call persp(y1(),s1())
    call persp(y2(),s2())
    plot lines: s1(1),s1(2); s2(1),s2(2)
    let x1(1) = ta
    let x1(2) = 0
    let x1(3) = tc
    let x2(1) = ta
    let x2(2) = 0
    let x2(3) = 0
    mat y1 = nn*x1
    mat y2 = nn*x2
    call persp(y1(),s1())
    call persp(y2(),s2())
    plot lines: s1(1),s1(2); s2(1),s2(2)
    let xx(1) = 0
    let xx(2) = 0
    let xx(3) = tc
    mat yy = nn*xx
    call persp(yy(),ss())
    

!   RAY OF SUNLIGHT AT NOON (ON THE DAY OF THE VERNAL EQUINOX)

    let delta = 10

    call polar (-delta,180,x3())
    mat x1 = eq*x3
    let tw = x1(3)/(x1(3)-tc)
    let x2(1) = (1-tw)*x1(1)
    let x2(2) = (1-tw)*x1(2)
    let x2(3) = 0
    mat y1 = nn*x1
    mat y2 = nn*x2
    call persp(y1(),s1())
    call persp(y2(),s2())
    set color 17    !   gold
    plot lines: s1(1),s1(2); s2(1),s2(2)
    box show sun1$ at s1(1)-2*rad,s1(2)-2*rad 
    box show zpole$ at ss(1)-rad,ss(2)-rad
    box show sun2$ at s2(1)-rad,s2(2)-rad

    let tr = 20

    call polar (-delta,180-tr,x3())
    mat x1 = eq*x3
    let tw = x1(3)/(x1(3)-tc)
    let x2(1) = (1-tw)*x1(1)
    let x2(2) = (1-tw)*x1(2)
    let x2(3) = 0
    mat y1 = nn*x1
    mat y2 = nn*x2
    call persp(y1(),s1())
    call persp(y2(),s2())
    set color 17    !   gold
    plot lines: s1(1),s1(2); s2(1),s2(2)

    call polar (-delta,180+tr,x3())
    mat x1 = eq*x3
    let tw = x1(3)/(x1(3)-tc)
    let x2(1) = (1-tw)*x1(1)
    let x2(2) = (1-tw)*x1(2)
    let x2(3) = 0
    mat y1 = nn*x1
    mat y2 = nn*x2
    call persp(y1(),s1())
    call persp(y2(),s2())
    set color 17    !   gold
    plot lines: s1(1),s1(2); s2(1),s2(2)

    for jj = -tr to tr
        call polar (-delta,180+jj,x3())
        mat x1 = eq*x3
        let tw = x1(3)/(x1(3)-tc)
        mat y1 = nn*x1
        call persp(y1(),s1())
        plot lines: s1(1),s1(2);
    next jj
    plot points: s1(1),s1(2)

    for jj = -tr to tr
        call polar (-delta,180+jj,x3())
        mat x1 = eq*x3
        let tw = x1(3)/(x1(3)-tc)
        let x2(1) = (1-tw)*x1(1)
        let x2(2) = (1-tw)*x1(2)
        let x2(3) = 0
        mat y2 = nn*x2
        call persp(y2(),s2())
        plot lines: s2(1),s2(2);
    next jj
    plot points: s2(1),s2(2)
    
    stop
    
!   SUBPROGRAMS

    sub persp(zz(),ss())

        let dp = eye/(eye - zz(3))
        let ss(1) = dp*zz(1)  
        let ss(2) = dp*zz(2)

    end sub

    sub unit(uu(),vv())

        call norm(uu(),a)
        for ii = 1 to 3
            let vv(ii) = uu(ii)/a
        next ii

    end sub

    sub norm(uu(),a)

        let b = dot(uu(),uu())
        let a = sqr(b)

    end sub

    sub cross(uu(),vv(),ww())

        let ww(1) = uu(2)*vv(3) - uu(3)*vv(2)
        let ww(2) = uu(3)*vv(1) - uu(1)*vv(3)
        let ww(3) = uu(1)*vv(2) - uu(2)*vv(1)

    end sub

    sub anti(zz(),cc(,))

        !  In practice, zz should be a unit vector.
        for ii = 1 to 3
            let cc(ii,ii) = 0
        next ii
        let cc(1,2) = -zz(3)
        let cc(2,1) =  zz(3)
        let cc(1,3) =  zz(2)
        let cc(3,1) = -zz(2)
        let cc(2,3) = -zz(1)
        let cc(3,2) =  zz(1)

    end sub

    sub exp(ff,tt,cc(,),ee(,))

        !  In practice, cc should be defined by a unit vector, so that cc^3=-cc.  See
        !  the foregoing subroutine.  Then ee = id + (1-cos(gg))*cc^2 + sin(gg)*cc =
        !  exp(gg*cc) where gg = 360*ff*tt.  One may prove that it is so by differentiating 
        !  with respect to gg.
        mat bb = cc*cc
        mat aa = id + bb
        let gg = 360*ff*tt
        let cosgg = cos(gg)
        mat dd = cosgg*bb
        mat bb = aa - dd
        let singg = sin(gg)
        mat dd = singg*cc
        mat ee = bb + dd
        
    end sub

    sub polar(th,ph,zz())

        let zz(1) = cos(th)*cos(ph)
        let zz(2) = cos(th)*sin(ph)
        let zz(3) = sin(th)

    end sub

end