#version 3.1

#include "colors.inc"

global_settings
{
  assumed_gamma 1.0
}


//====================================================
//  HERE'S THE CODE
//----------------------------------------------------
//
//  the major problem is that the parser chokes on the
//  macros.  Generally, I'v found that if you put a 
//  ";" after a macro call, the parser halts with some
//  kind of buggy error.  But if you don't use a ";",
//  it generates all kinds of error messages.  maybe
//  someone can clean up the code.
//====================================================


#declare maxstack = 1024;                       

//declare the first stack
#declare stackpos_1 = 0;
#declare stack_1 = array[maxstack][3]

//declare the second stack
#declare stackpos_2 = 0;
#declare stack_2 = array[maxstack][3]

//initialize them so the parser doesn't choke
#declare iterator = 0;
#while(iterator<maxstack)
        #declare stack_1[iterator][0] = <0,0,0>;
        #declare stack_1[iterator][1] = <0,0,0>;
        #declare stack_1[iterator][2] = <0,0,0>;
        #declare stack_2[iterator][0] = <0,0,0>;
        #declare stack_2[iterator][1] = <0,0,0>;
        #declare stack_2[iterator][2] = <0,0,0>;
        #declare iterator = iterator + 1;     
#end

//====================================================
//  define a bunch of macros to make my life easier.
//====================================================


#macro vector_string(vect)
        #local donothing = "<",str(vect.x,10,5),",",str(vect.y,10,5),",",str(vect.z,10,5),">";
        donothing
#end

//------------------------------------------------
//  STACK HANDLING ROUTINES
//------------------------------------------------

#macro push(s,sp,vbegin,vend,col)
        #if(sp+1<maxstack)
                #declare sp = sp + 1;
                #declare s[sp][0] = vbegin;
                #declare s[sp][1] = vend;
                #declare s[sp][2] = col;
                //#render concat("push :: vbegin ",vector_string(vbegin),"     vend :",vector_string(vend),"\n")
        #else
                #error concat("stack pointer :",str(sp,5,0),"  maxstack:",str(maxstack,5,0),"\n")
        #end
#end

#macro peek(s,sp,i)
        #if(sp>=0 & sp<maxstack)
                s[sp][i]
        #else
                #error "stack out of bounds"
        #end
#end

#macro pop(s,sp)
        #if(sp>0)
                #declare sp = sp - 1;
        #else
                #error "stack underflow"
        #end
#end


//------------------------------------------------
//  PERTURBATION ROUTINES
//------------------------------------------------

#macro perturb(vbegin,vend,s,amount)
        #local dv    = vend-vbegin;
        #local leng1 = vlength(dv);
        #local randv = <rand(s)-.5,rand(s)-.5,rand(s)-.5>*2*leng1*amount;
        #local leng2 = vlength(randv);  
        #local perpv = vnormalize(randv - vdot(randv,dv)/leng1/leng2 * dv)*leng1*amount; 
        #local vtemp = vend + perpv;
        vtemp
#end


#macro vmid(vbegin,vend,s,amount)
        #local vmiddle = vbegin + (vend-vbegin)*(0.5 + (rand(s)-0.5)*amount/2.0 );
        #local vtemp2 = perturb(vbegin,vmiddle,s,amount);
        vtemp2
        
        //#local vtemp = vbegin + (vpert-vbegin)*(0.5 + (rand(s)-0.5)*amount/2.0 );
        //#local leng  = vlength(vend-vbegin);
        //#local vtemp = vbegin + (vend-vbegin)*0.5 + <rand(s)-.5,rand(s)-.5,rand(s)-.5>*2*leng*amount;
#end


//------------------------------------------------
//  RECURSIVE SUBDIVISION ROUTINES
//------------------------------------------------


#macro split_stack(st_1,stp_1,st_2,stp_2,amount,fork,s)
        #declare mindepth = -1.0;
        #declare maxdepth = -1.0;
        #while (stp_1>0)                          
                #local vbegin  = peek(st_1,stp_1,0)
                #local vend    = peek(st_1,stp_1,1)
                #local cparam  = peek(st_1,stp_1,2)
                #local cdepth  = cparam.x;
                
                //ideally this is supposed to makr the lightning
                //dimmer after it forks, but it's not working right
                #local cintens = cparam.y;
                #local nintens = cintens/2;
                
                #if (mindepth < 0 )
                  #declare mindepth = cdepth.x;
                #end
                #if (maxdepth < 0 )
                  #declare maxdepth = cdepth.x;
                #end

                #declare mindepth = min(cdepth,mindepth);
                #declare maxdepth = max(cdepth,maxdepth);
                
                #if (cdepth<max_light)
                
                        #local vmiddle = vmid(vbegin,vend,s,amount)
                        
                        #local vleng_0 = vlength(vmiddle-vbegin);
                        #local depth_0 = -log(vleng_0/light_length)/log(2);
                        push(st_2,stp_2,vbegin,vmiddle,<depth_0,cintens,0>)
                        
                        #local vleng_1 = vlength(vend-vmiddle);
                        #local depth_1 = -log(vleng_1/light_length)/log(2);
                        
                        #if (rand(s)<fork*pow(2,-cdepth/2))
                                push(st_2,stp_2,vmiddle,vend,<depth_1,nintens,0>)
                                
                                #local length_fraction = vleng_0/(vleng_0 + vleng_1);
                                #local vend2 = perturb(light_begin,light_end,s,length_fraction*amount)
                                #local vleng_2 = vlength(vend2-vmiddle);
                                #local depth_2 = -log(vleng_2/light_length)/log(2);
                                push(st_2,stp_2,vmiddle,vend2,<depth_2,nintens,0>)
                        
                                #declare mindepth = min(min(min(depth_0,depth_1),depth_2),mindepth);
                                #declare maxdepth = max(max(max(depth_0,depth_1),depth_2),maxdepth);
                        #else
                                #declare mindepth = min(min(depth_0,depth_1),mindepth);
                                #declare maxdepth = max(max(depth_0,depth_1),maxdepth);

                                push(st_2,stp_2,vmiddle,vend,<depth_1,cintens,0>)
                        #end

                #else
                        push(st_2,stp_2,vbegin,vend,<cdepth,cintens,0>)
                #end
                
                #declare stp_1 = stp_1 - 1;
        #end
                
#end

//------------------------------------------------
//  DISPLAY ROUTINE
//------------------------------------------------


#macro display_stack(st,stp)
        #while(stp>0)
                #local v1    = peek(st,stp,0)
                #local v2    = peek(st,stp,1)
                #local param = peek(st,stp,2) 
                #local intens = param.y

                cylinder{ v1,v2,0.02  pigment {color rgb <1,1,1>*intens} }
                pop(st,stp)
        #end
#end


//------------------------------------------------
//  ACTUALLY DO THE CALCULATION
//------------------------------------------------

                                    
#declare light_begin = <0,1,0>;
#declare light_end = <0,0,0>;
#declare light_param = <0,1,0>;
#declare light_vect = light_end-light_begin;
#declare light_length = vlength(light_vect);

#declare light_fork = 0.25;
#declare light_dev  = 0.4;

#declare min_light = 3;
#declare max_light = 4;
#declare seed_1 = seed(-1);
#declare mindepth = -5.0;

push(stack_1,stackpos_1,light_begin,light_end,light_param)

split_stack(stack_1,stackpos_1,stack_2,stackpos_2,light_dev,light_fork,seed_1)
#render concat("***** stack pointer :",str(stackpos_2,5,0),"  mindepth:",str(mindepth,5,0),"\n")

#while(mindepth < min_light)
        split_stack(stack_2,stackpos_2,stack_1,stackpos_1,light_dev,light_fork,seed_1)
        #render concat("***** stack pointer :",str(stackpos_1,5,0),"  mindepth:",str(mindepth,5,0),"\n")
        split_stack(stack_1,stackpos_1,stack_2,stackpos_2,light_dev,light_fork,seed_1)
        #render concat("***** stack pointer :",str(stackpos_2,5,0),"  mindepth:",str(mindepth,5,0),"\n")
        
#end        


//------------------------------------------------
//  A VERY SIMPLE SCENE
//------------------------------------------------


camera
{
  location  <0.0, 0.5, -4.0>
  direction 1.5*z
  right     4/3*x
  look_at   <0.0, 0.0,  0.0>
}

sky_sphere
{
  pigment
  {
    gradient y
    color_map { [0.0 color blue 0.6] [1.0 color rgb 1] }
  }
}

light_source
{
  0*x // light's position (translated below)
  color red 1.0  green 1.0  blue 1.0  // light's color
  translate <-30, 30, -30>
}


display_stack(stack_2,stackpos_2)