rem Fractals 3D II ........ Rev 4.0
      rem A J Tooth // July 2007
      
      rem Preamble==================================
      on error if (err=17) then quit
      *FLOAT 64
      himem=lomem + 100000000
      install @lib$+"MyUtils.bbc"
      install @lib$+"BMP_Utils.bbc"
      rem Preamble==================================
      
      rem Setup
      proc_setup
      
      repeat
        rem Parameter reset
        proc_parset
        
        rem Initialise the transformations
        proc_setmat(ind&) : ind&=1
        
        for k&=0 to ml&-1
          
          rem Reset additive component
          proc_adj(k&)
          
          rem Reset colour
          Rr&=rnd(255) : Gg&=rnd(255) : Bb&=rnd(255)
          
          rem Draw next figure
          proc_figure
          
        next k&
        
        *REFRESH
        a$=inkey$(100) : cls
        
        origin 0,0
        rem Re display the plotted points
        for h&=1 to 150
          Dis+=2.0
          th+=0.5 : theta=pi*th/180
          ph+=0.5   : phi=pi*ph/180
          proc_redo(Dis,theta,phi)
        next h&
        
        *REFRESH
        a$=inkey$(100) : cls
      until a$<>""
      
      *REFRESH ON
      a$=get$
      quit
      rem End of Program ================================================
      rem ===============================================================
      
      rem Setup
      def proc_setup
      
      rem Maximise display
      proc_fullscreen(xscreen%,yscreen%)
      
      Num%=70000 : ml&=3 : Lim%=ml&*Num%
      dim poynt{(ml&*Num%)vert{x,z,y},cll{ro&,go&,bo&,r&,g&,b&}}
      
      rem Set up a screen bitmap
      proc_BMP_Set(xscreen%,yscreen%,pic%,Wlim%,lgth%)
      
      dim Pltt% 1000, Clear% 200, Poynt% 300, colset% 500, itmp% 3, ftmp% 7
      dim calc% 500, Master% 200
      
      rem Dual-pass assembly, in case of labels
      for pass&=0 to 2 step 2
        proc_Master(pass&)
        proc_Pltt(pass&)
        proc_Clear(pass&)
        proc_Poynt(pass&)
        proc_colset(pass&)
        proc_calc(pass&)
      next pass&
      
      *REFRESH OFF
      
      ind&=0 : Flp&=0
      endproc
      rem ===============================================================
      
      rem Parameter reset
      def proc_parset
      
      A=-2.0 : rem Camera position
      Dis=1000.0 : rem Base distance
      th=-0.5 : ph=-0.5 : rem Initial camera angles
      Cnt%=0 : rem Count the number of points plotted
      origin xscreen%,yscreen%
      endproc
      rem ===============================================================
      
      rem Draw next figure
      def proc_figure
      local i%,j&,xnew,ynew,znew,p&,pic
      
      for i%=1 to Num%
        Cnt%+=1
        pic=rnd(1)
        if pic<=prob(0) then p&=0
        if pic>prob(0) and pic<prob(1) then p&=1
        if pic>prob(1) and pic<prob(2) then p&=2
        if pic>prob(2) then p&=3
        
        xnew=xy(0)*mat(0,0,p&) + xy(1)*mat(0,1,p&) + xy(2)*mat(0,2,p&) + adj(0,p&)
        znew=xy(0)*mat(1,0,p&) + xy(1)*mat(1,1,p&) + xy(2)*mat(1,2,p&) + adj(1,p&)
        ynew=xy(0)*mat(2,0,p&) + xy(1)*mat(2,1,p&) + xy(2)*mat(2,2,p&) + adj(2,p&)
        
        xy(0)=xnew : xy(1)=znew : xy(2)=ynew
        
        rem Convert from 3D to 2D screen coords
        proc_poynt(Dis,xy(0),xy(1),xy(2),Px%,Py%)
        
        poynt{(Cnt%)}.vert.x=xy(0)
        poynt{(Cnt%)}.vert.z=xy(1)
        poynt{(Cnt%)}.vert.y=xy(2)
        
        poynt{(Cnt%)}.cll.ro&=Rr&
        poynt{(Cnt%)}.cll.go&=Gg&
        poynt{(Cnt%)}.cll.bo&=Bb&
        
        rem Alter the hue
        proc_hue(1,Cnt%,xy(1),poynt{})
        
        move Px%,Py% : draw Px%,Py%
        
        if i%mod500=0 then *REFRESH
        
      next i%
      endproc
      rem ===============================================================
      
      rem Alter the hue
      def proc_hue(Ind&,i%,Znn,poynt{})
      local red,R&,G&,B&
      
      rem Set the current hue
      call colset%
      
      if Ind&=1 then colour 2,R&,G&,B& : gcol 2
      endproc
      rem ===============================================================
      
      rem Set transformation matrix
      def proc_setmat(Ctr&)
      local i&,Fac,Theta,a&,b&
      
      if Ctr&=0 then
        dim mat(2,2,3), xy(2), prob(3)
      endif
      
      xy(0)=500.0
      xy(1)=0.0
      xy(2)=500.0
      
      for i&=0 to 3
        case Flp& of
          when 0:
            for a&=0 to 2
              for b&=0 to 2
                mat(a&,b&,i&)=(2*rnd(1)-1)/1.5
              next b&
            next a&
          when 1:
            Fac=rnd(1)
            Th=2*pi*rnd(1) : Ph=pi*rnd(1)-pi/2 : Ps=2*pi*rnd(1)
            
            mat(0,0,i&)=Fac*( cos(Ps)*cos(Ph)*cos(Th) - sin(Ps)*sin(Th))
            mat(0,1,i&)=Fac*(-sin(Ps)*cos(Ph)*cos(Th) - cos(Ps)*sin(Th))
            mat(0,2,i&)=Fac*(-sin(Ph)*cos(Th))
            mat(1,0,i&)=Fac*( cos(Ps)*cos(Ph)*sin(Th) + sin(Ps)*cos(Th))
            mat(1,1,i&)=Fac*(-sin(Ps)*cos(Ph)*sin(Th) + cos(Ps)*cos(Th))
            mat(1,2,i&)=Fac*(-sin(Ph)*sin(Th))
            mat(2,0,i&)=Fac*( cos(Ps)*sin(Ph))
            mat(2,1,i&)=Fac*(-sin(Ph)*sin(Ps))
            mat(2,2,i&)=Fac*( cos(Ph))
            
        endcase
      next i&
      
      prob(0)=0.5*rnd(1)
      prob(1)=prob(0)+0.25*rnd(1)
      prob(2)=prob(1)+0.25*rnd(1)
      prob(3)=1.0
      
      Flp&=1-Flp&
      endproc
      rem ===============================================================
      
      rem Reset additive component
      def proc_adj(Ctr%)
      local i&
      
      if Ctr%=0 then
        dim adj(2,3)
      endif
      
      for i&=0 to 3
        adj(0,i&)=1.0*(rnd(500) - 250)
        adj(1,i&)=1.0*(rnd(500) - 250)
        adj(2,i&)=1.0*(rnd(500) - 250)
      next i&
      endproc
      rem ===============================================================
      
      def proc_poynt(Ds,Xx,Zz,Yy,return Px%,return Py%)
      local fac
      
      fac=1.0/(1.0-(Ds+Zz)/A) : rem Reduction or scale factor
      Px%=int(xscreen%*Xx*fac) : Py%=int(yscreen%*Yy*fac)
      
      endproc
      rem ===================================================================
      
      rem Re display the plotted points
      def proc_redo(D,th%,ph%)
      local i%,X,Z,Y
      
      call Master%
      
      proc_BMP_DispB(xscreen%,yscreen%,pic%)
      *REFRESH
      rem Clear the BMP to zeros
      call Clear%
      
      endproc
      rem =======================================================================
      
      rem Master Routine
      def proc_Master(opt&)
      P%=Master%
      [opt opt&
      
      mov edx,1
      mov [^i%],edx
      
      .iloop
      call calc%
      
      call Poynt%
      
      call colset%
      
      call Pltt%
      
      inc dword [^i%]
      mov edx,[^i%]
      cmp edx,[^Lim%]
      jbe near iloop
      
      ret
      ]
      endproc
      rem ======================================================================
      
      rem BMP Update Routine
      def proc_Pltt(opt&)
      P%=Pltt%
      [opt opt&
      
      mov eax,[^i%]
      imul eax,30
      mov cl,(^poynt{(0)}.cll.r&)[eax]
      mov [^Rn&],cl
      mov cl,(^poynt{(0)}.cll.g&)[eax]
      mov [^Gn&],cl
      mov cl,(^poynt{(0)}.cll.b&)[eax]
      mov [^Bn&],cl
      
      mov ebx,[^Py%]
      add ebx,[^yscreen%]
      shr ebx,1
      cmp ebx,0
      jl near miss
      cmp ebx,[^yscreen%]
      jge near miss
      
      mov ebx,[^Px%]
      add ebx,[^xscreen%]
      shr ebx,1
      cmp ebx,0
      jl near miss
      cmp ebx,[^xscreen%]
      jge near miss
      
      mov eax,[^Py%]
      add eax,[^yscreen%]
      shr eax,1
      imul eax,[^Wlim%]
      mov ebx,[^Px%]
      add ebx,[^xscreen%]
      shr ebx,1
      add eax,ebx
      add eax,ebx
      add eax,ebx
      add eax,54
      
      mov cl,[^Bn&]
      mov pic%[eax],cl
      mov cl,[^Gn&]
      mov pic%[eax+1],cl
      mov cl,[^Rn&]
      mov pic%[eax+2],cl
      
      .miss
      ret
      ]
      endproc
      rem ======================================================================
      
      rem BMP Update Routine
      def proc_Clear(opt&)
      P%=Clear%
      [opt opt&
      
      mov edx,54
      mov [^cc%],edx
      
      mov cl,0
      .ccloop
      mov eax,[^cc%]
      mov pic%[eax],cl
      
      inc dword [^cc%]
      mov edx,[^cc%]
      cmp edx,[^lgth%]
      jb near ccloop
      
      ret
      ]
      endproc
      rem ======================================================================
      
      rem 3D to 2D conversion
      def proc_Poynt(opt&)
      P%=Poynt%
      [opt opt&
      
      finit
      fld qword [^D]            ;fac=1.0/(1.0-(D+Znn)/A)
      fld qword [^Znn]
      faddp st1,st0
      fld qword [^A]
      fdivp st1,st0
      fld1
      fsubrp st1,st0
      fld1
      fdivrp st1,st0            ;Calc 0f fac
      fld st0
      
      fild dword [^xscreen%]
      fld qword [^Xnn]
      fmulp st1,st0
      fmulp st1,st0
      fistp dword [^Px%]
      
      fild dword [^yscreen%]
      fld qword [^Ynn]
      fmulp st1,st0
      fmulp st1,st0
      fistp dword [^Py%]
      
      ret
      ]
      endproc
      rem ======================================================================
      
      rem Colour set
      def proc_colset(opt&)
      P%=colset%
      [opt opt&
      
      finit
      fild dword [i500%]
      fld qword [^Znn]
      fsubp st1,st0
      fild dword [i500%]
      fdivp st1,st0
      fstp qword [^red]
      
      mov eax,[^i%]
      imul eax,30
      mov cl,(^poynt{(0)}.cll.ro&)[eax]
      mov [^RO&],cl
      mov cl,(^poynt{(0)}.cll.go&)[eax]
      mov [^GO&],cl
      mov cl,(^poynt{(0)}.cll.bo&)[eax]
      mov [^BO&],cl
      
      mov ebx,0
      mov bl,[^RO&]
      mov [itmp%],ebx
      fild dword [itmp%]
      fld qword [^red]
      fmulp st1,st0
      fistp dword [itmp%]
      
      mov ebx,[itmp%]
      cmp ebx,255
      jle okr1
      mov bl,255
      mov [^R&],bl
      jmp nxt1
      .okr1
      cmp ebx,0
      jge okr2
      mov bl,0
      mov [^R&],bl
      jmp nxt1
      .okr2
      mov [^R&],bl
      
      .nxt1
      mov ebx,0
      mov bl,[^GO&]
      mov [itmp%],ebx
      fild dword [itmp%]
      fld qword [^red]
      fmulp st1,st0
      fistp dword [itmp%]
      
      mov ebx,[itmp%]
      cmp ebx,255
      jle okg1
      mov bl,255
      mov [^G&],bl
      jmp nxt2
      .okg1
      cmp ebx,0
      jge okg2
      mov bl,0
      mov [^G&],bl
      jmp nxt2
      .okg2
      mov [^G&],bl
      
      .nxt2
      mov ebx,0
      mov bl,[^BO&]
      mov [itmp%],ebx
      fild dword [itmp%]
      fld qword [^red]
      fmulp st1,st0
      fistp dword [itmp%]
      
      mov ebx,[itmp%]
      cmp ebx,255
      jle okb1
      mov bl,255
      mov [^B&],bl
      jmp nxt3
      .okb1
      cmp ebx,0
      jge okb2
      mov bl,0
      mov [^B&],bl
      jmp nxt3
      .okb2
      mov [^B&],bl
      
      .nxt3
      mov eax,[^i%]
      imul eax,30
      mov cl,[^R&]
      mov (^poynt{(0)}.cll.r&)[eax],cl
      mov cl,[^G&]
      mov (^poynt{(0)}.cll.g&)[eax],cl
      mov cl,[^B&]
      mov (^poynt{(0)}.cll.b&)[eax],cl
      
      jmp ffin
      .i500% DD 500
      .ffin ret
      ]
      endproc
      rem ======================================================================
      
      rem Rotation Calculation
      def proc_calc(opt&)
      P%=calc%
      [opt opt&
      
      mov eax,[^i%]
      imul eax,30
      
      finit
      fld qword [^theta]
      fsincos
      fld qword (^poynt{(0)}.vert.x)[eax]
      fmulp st1,st0
      fld qword (^poynt{(0)}.vert.z)[eax]
      fmulp st2,st0
      fsubrp st1,st0
      fstp qword [^Xn]
      
      fld qword [^theta]
      fsincos
      fld qword (^poynt{(0)}.vert.z)[eax]
      fmulp st1,st0
      fld qword (^poynt{(0)}.vert.x)[eax]
      fmulp st2,st0
      faddp st1,st0
      fstp qword [^Zn]
      fld qword (^poynt{(0)}.vert.y)[eax]
      fstp qword [^Yn]
      
      fld qword [^phi]
      fsincos
      fld qword [^Zn]
      fmulp st1,st0
      fld qword [^Yn]
      fmulp st2,st0
      fsubrp st1,st0
      fstp qword [^Znn]
      
      fld qword [^phi]
      fsincos
      fld qword [^Yn]
      fmulp st1,st0
      fld qword [^Zn]
      fmulp st2,st0
      faddp st1,st0
      fstp qword [^Ynn]
      fld qword [^Xn]
      fstp qword [^Xnn]
      
      ret
      ]
      endproc
      rem ======================================================================
      
      def proc_dump
      
      Xn=X*cos(theta) - Z*sin(theta)
      Zn=X*sin(theta) + Z*cos(theta)
      Yn=Y
      
      Xnn=Xn
      Znn=Zn*cos(phi) - Yn*sin(phi)
      Ynn=Zn*sin(phi) + Yn*cos(phi)
      endproc
      rem ======================================================================