// TWIGS3D / Xmas Trees '07.11
// by Tsutomu HIGO  URL:http://www.asahi-net.or.jp/~nj2t-hg/

#version 3.6;
global_settings {
  assumed_gamma 2.2
  }

light_source {<50, 200, -100> color <1, 1, 1>}

camera {
  location <-15, 5, -40>
  angle 100
  look_at <10, 6, 20>
  }

#declare R1 = seed(3);

#macro Dla36 (Nmax, Rds, Ar)
union {
  #declare Cp = array[Nmax][7]
  #declare Cp[0][0] = 0; // x of sphere's center
  #declare Cp[0][1] = 0; // z of sphere's center
  #declare Cp[0][2] = 0; // y of sphere's center
  #declare Cp[0][3] = 0; // Distanse of sphere's center
  #declare Cp[0][4] = 2*pi*rand(R1); // Longitude Angle of sphere's center (keido)
  #declare Cp[0][5] = pi/2*(1-Ar+Ar*rand(R1)); // Latitude Angle of sphere's center (ido)
  #declare Cp[0][6] = 1; // Terminal Flag
  sphere {
    <Cp[0][0], Cp[0][2], Cp[0][1]> Rds
    pigment {color rgb <1, 1-0/Nmax*0.8, 0.2>}
    finish {
      reflection 0.3
      phong 2
      brilliance 1
      ambient 0.4
      }
    no_shadow
    }

  #declare Np = 1; #while (Np < Nmax)
    #declare Cp[Np][4] = 2*pi*rand(R1);
    #declare Cp[Np][5] = pi/2*(1-Ar+Ar*rand(R1));
    #declare Cp[Np][6] = 1;
    #if (Np = 1)
          #declare Cp[Np][3] = 2;
          #declare Cp[Np][0] = Cp[Np][3]*cos(Cp[Np][5])*cos(Cp[Np][4]);
          #declare Cp[Np][1] = Cp[Np][3]*cos(Cp[Np][5])*sin(Cp[Np][4]);
          #declare Cp[Np][2] = Cp[Np][3]*sin(Cp[Np][5]);
          #declare Nsclose = 0;
    #else
      #declare Ns = Np-1; #while (Ns>0)
          #declare Cp[Np][3] = Cp[Ns][3]; //Templary size
          #declare Cp[Np][0] = Cp[Np][3]*cos(Cp[Np][5])*cos(Cp[Np][4]);
          #declare Cp[Np][1] = Cp[Np][3]*cos(Cp[Np][5])*sin(Cp[Np][4]);
          #declare Cp[Np][2] = Cp[Np][3]*sin(Cp[Np][5]);
          #declare Det = 0;
          #declare Ldis = sqrt(pow(Cp[Np][0]-Cp[Ns][0],2)+pow(Cp[Np][1]-Cp[Ns][1],2)+pow(Cp[Np][2]-Cp[Ns][2],2));
        #if (Ldis < 1)
          #declare Apn = asin(Ldis/2/Cp[Ns][3]);
          #declare Cp[Np][3] = Cp[Ns][3]*cos(2*Apn)+sqrt(4-pow(Ldis*cos(Apn),2)); // Real size
          #declare Cp[Np][0] = Cp[Np][3]*cos(Cp[Np][5])*cos(Cp[Np][4]);
          #declare Cp[Np][1] = Cp[Np][3]*cos(Cp[Np][5])*sin(Cp[Np][4]);
          #declare Cp[Np][2] = Cp[Np][3]*sin(Cp[Np][5]);
          #declare Nsclose = Ns;
          #declare Det = 1;
          #declare Ns = 1;
        #end
      #declare Ns = Ns-1; #end
        #if (Det = 0)
          #declare Cp[Np][3] = 2;
          #declare Cp[Np][0] = Cp[Np][3]*cos(Cp[Np][5])*cos(Cp[Np][4]);
          #declare Cp[Np][1] = Cp[Np][3]*cos(Cp[Np][5])*sin(Cp[Np][4]);
          #declare Cp[Np][2] = Cp[Np][3]*sin(Cp[Np][5]);
          #declare Nsclose = 0;
        #end
    #end
          #declare Cp[Nsclose][6] = 0;
    union {
      sphere {<Cp[Np][0], Cp[Np][2], Cp[Np][1]> Rds}
      cylinder {<Cp[Np][0], Cp[Np][2], Cp[Np][1]> <Cp[Nsclose][0], Cp[Nsclose][2], Cp[Nsclose][1]> Rds}
      pigment {color rgb <0.7-Np/Nmax*0.7, 0.7, 0.3-Np/Nmax*0.2>}
      finish {
        reflection 0.3
        phong 2
        brilliance 1
        ambient 0.4
        }
      no_shadow
      }
  #declare Np = Np+1; #end

  #declare Np = 1; #while (Np < Nmax) 
    #if (Cp[Np][6] = 1)
      sphere {
        <Cp[Np][0], Cp[Np][2], Cp[Np][1]> Rds*3
        pigment {color rgbt <1, 1, 1, 0.95>}
        finish {
          reflection 0.7
          phong 2
          brilliance 1
          ambient 0.4
          }
        }
    #end
  #declare Np = Np+1; #end

  rotate y*-130
  }

#end

Dla36 (1400, 0.1, 1)

#declare Nt = -3; #while (Nt < 8)
  union {Dla36 (100, 0.1, 0.4) translate <30, 0, 0> rotate y*-30*Nt}
#declare Nt = Nt+1; #end

plane {
  y, -1.0
  hollow on
  pigment {color rgbt <1, 1, 1,0.8>}
  finish {
    reflection 1
    phong 2
    brilliance 1
    ambient 0.4
    }
  }

background {color rgb <0.01, 0, 0.2>}

fog {
  fog_type 2
  distance 300
  color rgb <0.4, 0.3, 0.8>
  fog_offset 70
  fog_alt 0.2
  turbulence 0.6
  }