Hi,

I'm trying to use the isosurface functionality of POV-Ray to
depict (in a loose, not physically accurate, way ;) an isosurface
of magnetic field strength around a moving point charge.

I'm using the "Biot-Savart law" where, if the particle has
velocity v and you put the particle at the origin, the field
strength at some point p is proportional to | ((v x p) / |p|^3) |
(where 'x' is the vector cross-product, and '|z|' means "vector
length of z"). 

I've worked around the difficulty of not being able to access
vector functions inside POV-Ray function definitions by defining
vector length and "cross product magnitude" functions, but I still
can't get it to work. Instead of a nice donut facing in the
direction of v, I just get whatever bounding box I choose for
the isosurface, completely filled up.

Here's my povray source

8<---C-U-T--H-E-R-E----
#version 3.5;
#include "colors.inc"
#include "functions.inc"
#declare CamLoc = < -28000000, 0,0>;
camera { location CamLoc look_at <0,0,0> angle 45 }
light_source {CamLoc color White*0.35}
light_source {<-50, 150,-75> color White}
#declare IsoFinish = finish { ambient 1 diffuse 1 specular 1
                              roughness 0.02 brilliance 2 }
#declare vel = <35,25,45>;
#declare velx = vel.x;
#declare vely = vel.y;
#declare velz = vel.z;
#declare f_vlen = function(x,y,z) { sqrt(x^2 + y^2 + z^2) }
#declare f_vcrossmag = function(v1x,v1y,v1z,v2x,v2y,v2z) {
        sqrt( (v1y*v2z - v1z*v2y)^2 + (v1z*v2x - v1x*v2z)^2
              + (v1x*v2y - v1y*v2x)^2 )
}
#declare f_bfield = function(x,y,z,vx,vy,vz) {
        f_vcrossmag(vx,vy,vz,x,y,z) / pow(f_vlen(x,y,z),3)
}
isosurface { function {f_bfield(x,y,z,velx,vely,velz) - 6}
        max_gradient 35000 threshold 0 contained_by {box{-8000000, 8000000}}
        texture { pigment {color Red} finish {IsoFinish} }
        scale 1/vlength(1)
}
8<---C-U-T--H-E-R-E----

Being convinced that the equation should yield more sensible
results, I wrote a perl script that positions spheres
along vectors chosen at random, attempting to find the point
along the vector for which the equation gives a magnitude
of 6 for the field, so it does a very rough approximation of
what I expected the POV-Ray isosurface to produce. In order
to make this a good debugging test, I implemented the functions
in perl in exactly the same way I've done them in povray.

This seems to work as expected, giving the plump donut shape,
lined up with the velocity vector.

As I'm by no means fully conversant with either the maths or the
POV-Ray isosurfaces and functions, it's quite possible I've
made a mistake somewhere - perhaps in my basic understanding of
what isosurfaces do and their limits, which is why I'm posting
this in the newusers group - but I *have* checked over the
equivalence between my perl and POV-Ray definitions of the
functions quite carefully, so I'd be interested if anyone can
tell me where my POV-Ray script is going wrong!

Just for reference, here's the perl script I used - if you want
to run it, you'll need to send the output to a file, like
this:
   ./bf.pl > bfield.pov

8<----C-U-T--H-E-R-E----
#! /usr/bin/perl
$velx = 35; $vely = 25; $velz = 45; $cmag = 6; $points = 10000;
print STDERR "velocity <$velx,$vely,$velz>, magnitude $cmag.\n";
sub vlen($$$) {
        my($x,$y,$z) = @_;
        sqrt($x*$x + $y*$y + $z*$z);
}
sub vcrossmag($$$$$$) {
        my ($ux,$uy,$uz,$vx,$vy,$vz) = @_;
        sqrt(($uy*$vz - $uz*$vy)**2 + ($uz*$vx - $ux*$vz)**2
             + ($ux*$vy - $uy*$vx)**2);
}
sub f_bfield($$$$$$) {
        my ($x,$y,$z,$vx,$vy,$vz) = @_; 
        vcrossmag($vx,$vy,$vz,$x,$y,$z) / ((vlen($x,$y,$z))**3);
}
sub findmag($$$$) {
        my ($x,$y,$z,$mag) = @_;
        my $scale = 1; my $tried = 0; my $bestdiff = $mag;
        my $bestscale = $scale;
        while($tried++ <= 1000) {
                $m = f_bfield($x*$scale,$y*$scale,$z*$scale,$velx,$vely,$velz);
                if (abs($m - $mag) < $bestdiff) {
                        $bestdiff = abs($m - $mag);
                        $bestscale = $scale;
                }
                if ($m == 0) {return 0;}
                if (abs($m-$mag) <= $mag/100) {return $scale;}
                $scale *= sqrt($m/$mag);
        }
        print STDERR "kludge not working, returning $bestscale with bestdiff
$bestdiff\n";
        return $bestscale;
}
$vvscale = vlen($velx,$vely,$velz);
$vvx = $velx/$vvscale;
$vvy = $vely/$vvscale;
$vvz = $velz/$vvscale;
print "#declare velvect = <$vvx,$vvy,$vvz>;\n";
print "#declare velscale = $vvscale;\n";
print "#declare bobbins = union {\n";
for (1..$points) {
        $x = rand(0.2) - 0.1;
        $y = rand(0.2) - 0.1;
        $z = rand(0.2) - 0.1;
        $scale = findmag($x,$y,$z,$cmag);
        ($x,$y,$z) = (($x*$scale), ($y*$scale), ($z*$scale));
        print "\tsphere { <$x,$y,$z>, 0.05 }\n";
}
print "}\n";
print <<END_POV;
#include "colors.inc"
global_settings {
        assumed_gamma 2.2
}
cylinder { 0, velvect*velscale 0.1 pigment{rgb<1,1,1>}}
#declare CamLoc = vnormalize(vcross(y,vcross(x,velvect))) * velscale/5;
camera {
        location CamLoc
        look_at < 0, 0, 0>
//        direction <0,0,1.3333>
        angle 45
}
light_source {CamLoc color White*0.35}
light_source {<-50, 150,-75> color White}
#declare IsoFinish =
finish {
        ambient 0.31 diffuse 1
        specular 1 roughness 0.02
        brilliance 2
}
object {
        bobbins
        rotate y*clock*180
        pigment { rgb<0.9,0.1,0.1>}
        finish {IsoFinish}
}
END_POV
8<----C-U-T--H-E-R-E----

Well, I hope that's not been *too* long a post, and thanks for having
a look!

- Pete.

Post a reply to this message