"Simon de Vet"  wrote
> This is what I want POV to do. I think I can construct the Nth zone
without
> much difficulty if I know the locations of the Nth nearest points, but I
> don't know how hard this will be to do. Any ideas?

I wrote the code for Nth Brillouin zones in Fortran (and later in C)
long time (15 years?) ago.  The original purpose of the development
of the isosurface patch was to visualize the Brillouin zones and
Fermi surfaces (=isosurfaces in BZ).
However, it might be difficult to use those data with "isosurface"
in POV-Ray 3.5 because of the lack of "data_*D_*" functions.

Anyway, you can visualize the zones with only POV code.
The below code is an example of Brillouin zones of
fcc crystals. (it's messy, though)

R. Suzuki


// The code
#version 3.5;

#include "colors.inc"
#include "functions.inc"
#include "woods.inc"

#declare N_BZ=1;       // n-th BZ     1st - 4th
#declare N_POINTS=200;  // higher value : sharper edge but slower

#declare BT=array[130][3]
#declare AL=array[130]
#declare  BZ=array[12]

#declare L=0;
#declare i=0;
#while(i<3)
  #declare j=0;
  #while(j<3)
    #declare k=0;
    #while(k<3)
         #declare ZZ=0;
         #if (int(i)=1)
          #if (int(j)=1)
           #if (int(k)=1)
             #declare ZZ=1;
             #debug "zero"
           #end
          #end
         #end

         #if (ZZ=0)
              #declare BT[L][0]=-(i-1);
              #declare BT[L][1]=-(j-1);
              #declare BT[L][2]=-(k-1);
              #declare
AL[L]=sqrt(BT[L][0]*BT[L][0]+BT[L][1]*BT[L][1]+BT[L][2]*BT[L][2])*0.5;
              #debug concat(str(AL[L],6,2),str(BT[L][2],6,2),"\n")
              #declare BT[L][0]=BT[L][0]/AL[L]*0.5;
              #declare BT[L][1]=BT[L][1]/AL[L]*0.5;
              #declare BT[L][2]=BT[L][2]/AL[L]*0.5;
              #declare L=L+1;
         #end
         #declare k=k+1;
      #end
      #declare j=j+1;
    #end
    #declare i=i+1;
  #end
  #debug "\n"

#declare i=0;
#while(i<2)
  #declare j=0;
  #while(j<2)
    #declare k=0;
    #while(k<2)
              #declare BT[L][0]=-(i-0.5);
              #declare BT[L][1]=-(j-0.5);
              #declare BT[L][2]=-(k-0.5);
              #declare
AL[L]=sqrt(BT[L][0]*BT[L][0]+BT[L][1]*BT[L][1]+BT[L][2]*BT[L][2])*0.5;
              #debug concat(str(AL[L],6,2),str(BT[L][2],6,2),"\n")
              #declare BT[L][0]=BT[L][0]/AL[L]*0.5;
              #declare BT[L][1]=BT[L][1]/AL[L]*0.5;
              #declare BT[L][2]=BT[L][2]/AL[L]*0.5;
              #declare L=L+1;
         #declare k=k+1;
      #end
      #declare j=j+1;
    #end
    #declare i=i+1;
  #end
  #debug "\n"


#declare BB=object {blob{
 threshold 1

#declare i=-0.01;
#while (i<(pi*0.5+0.001))
  #declare j=pi*0.25;
  #while (j<(pi*0.5+0.001))
      #declare XX=sin(i)*cos(j);
      #declare YY=sin(i)*sin(j);
      #declare k=0;
      #declare BZ[0]=1e10;
      #declare BZ[1]=1e10;
      #declare BZ[2]=1e10;
      #declare BZ[3]=1e10;
      #while (k<L)
         #declare COSAL= BT[k][0]*XX+BT[k][1]*YY+BT[k][2]*cos(i);
         #declare LENGTH=AL[k]/(COSAL+0.00000001);
         #if (LENGTH>0)
           #if (BZ[0]>LENGTH)
            #declare BZ[3]=BZ[2];
            #declare BZ[2]=BZ[1];
            #declare BZ[1]=BZ[0];
            #declare BZ[0]=LENGTH;
           #else
             #if (BZ[1]>LENGTH)
               #declare BZ[3]=BZ[2];
               #declare BZ[2]=BZ[1];
               #declare BZ[1]=LENGTH;
             #else
               #if (BZ[2]>LENGTH)
                 #declare BZ[3]=BZ[2];
                 #declare BZ[2]=LENGTH;
               #else
                 #if (BZ[3]>LENGTH)
                   #declare BZ[3]=LENGTH;
                 #end
               #end
             #end
           #end
         #end
         #declare k=k+1;
      #end
      #declare COSAL=BZ[N_BZ-1];
      sphere{<XX*COSAL*2, YY*COSAL*2, cos(i)*COSAL*2> 12/N_POINTS, 1}
      #if (j>pi*0.25)
        sphere{<YY*COSAL*2, XX*COSAL*2, cos(i)*COSAL*2> 12/N_POINTS, 1}
      #end
      #declare j=j+pi/N_POINTS;
    #end
    #declare i=i+pi/N_POINTS;
     #debug concat(str(i,6,3),"\n")
  #end

}}


// ----------------------------------------
global_settings {
  assumed_gamma 1.0
}


camera {
  location  <0.0, 1.5, -4.0>
  direction 2.5*z
  right     x*image_width/image_height
  look_at   <0.0, 0.5,  0.0>
}

sky_sphere {
  pigment {
    gradient y
    color_map {
      [0.0 rgb <0.4,0.5,0.35>]
      [0.7 rgb <0.02,0.05,0.1>]
    }
  }
}

light_source {
  <0, 0, 0>            // light's position (translated below)
  color rgb <1, 1, 1>  // light's color
  translate <-30, 30, -30>
}


#declare CC=object{
union{
  object {BB}
  object {BB rotate z*90}
  object {BB rotate z*180}
  object {BB rotate z*270}
}
}

object{
  union{
    object {CC}
    object {CC scale <1,1,-1> }
  }
  texture{T_Wood4 scale 0.5}
  scale 0.4
  translate y*0.5
  rotate y*-20
}

// end of the code

Post a reply to this message