Print

Please login or register.
1 Hour 1 Day 1 Week 1 Month Forever
Login with username, password and session length

[da_nang]

Hi, who just loves conics today? Even better, solid angles!
Yeah, the Monte Carlo approach, while seemingly accurate to a degree, was seemingly too slow for my purposes. So I tried numerical integration instead. And this function is putting out nonsense that is either contributed to programmer error or numerical instability. I can't tell which.

Code: [Select]

function [p] = generateProbSeq(c, b, r, d, t)
%GENERATEPROBSEQ Generates probability sequence for initial CoF c,
%   bloom b and max rate of fire r at range d for fire pattern t.
p = zeros(length(t),1);
alpha = c*pi/180;
dt = [diff(t);0];
dir = [0 0 -1]';



    function [prob] = fastProb(d, alpha, dir)
        %Computes the probability by calculating the solid angle of the part
        %of the rectangle inside the cone.
        % d is the distance from the xy plane, alpha is the apex angle and
        % dir is the direction of the cone. The cone must be above the xy
        % plane.
       
        %Main integrand
        fun = @(x,y) d./(x.^2 + y.^2 + d.^2).^(1.5);
       
        if d < 0
            %Apex below xy plane
            prob = 0;
            return;
        elseif d == 0
            %Point blank
            phi = acos(dot([0 0 1],dir)/sqrt(sum(dir.^2)));%azimuth
            if phi > pi/2 - alpha
                %Part of the cone is still intersecting
                prob = 1;
            else
                %No intersection.
                prob = 0;
            end
            return;
        end
       
        len = 1.8;%m
        wid = 0.88/1.8;%m
       
        h = [0 0 d]';
       
        z1 = dir/sqrt(sum(dir.^2));%Normalize
        A = null(z1');
        x1 = A(:,1);%Orthonormal to z
        y1 = A(:,2);%Orthonormal to z
        %Both orthonormal to each other
       
        if alpha == 0 %Line must intersect with rectangle for hit
            k_z = (-h'*[0 0 1]')/(z1'*[0 0 1]');
            point = h + k_z*z1;
            if k_z > 0 && abs(point(1,1)) <= len/2 && abs(point(2,1)) <= wid/2
                %Is pointing towards it and intersects rectangle
                prob = 1;
            else
                %Is not
                prob = 0;
            end
            return;
        end
       
        %Solid angle of spherical cap
        solidAngleCone = 2*pi*(1-cos(alpha/2));
       
        %Solid angle of rectangle
        solidAngleRect = dblquad(@(x,y) ...
            isInCone(x,y,h, x1, y1, z1, alpha).*fun(x,y),...
            -len/2, len/2, -wid/2, wid/2);
       
        %Final probability
        %Should be between 0 and 1
        prob = solidAngleRect/solidAngleCone;
       
               
    end


for i = 1:length(t);
    [prob] = fastProb(d, alpha, dir);
    p(i) = prob;   
    alpha = alpha + b*pi/180;
    if (dt(i) > 60/r + 0.05)
        alpha = c*pi/180;
    end
end
end

function [yes] = isInCone(x, y, h, x1, y1, z1, alpha)
        %x is a vector, y a scalar
        [m, ~] = size(x);
        if (m ~= 1)
            xs = x';
        else
            xs = x;
        end
       
        A = zeros(3, length(xs));
        A(1,:) = xs;
        A(2,:) = y*ones(1,length(xs));
        X = repmat(x1, 1, length(xs));
        Y = repmat(y1, 1, length(xs));
        Z = repmat(z1, 1, length(xs));
        H = repmat(h, 1, length(xs));
       
        x2 = dot(A-H,X);
        y2 = dot(A-H,Y);
        z2 = dot(A-H,Z);
       
        yes = (x2.^2 + y2.^2 <= ...
            ((z2.^2).*tan(alpha/2)).^2) & (z2 >= 0);

        if isequal(size(yes), size(x)) == 0
            yes = yes';
        end
end




[MagmaMcFry]

Okay, I must admit I don't really have the patience to walk through your opaque code. But I would like to suggest a solution that runs in O(1) time and needs no approximations or integration at all, assuming you're still working on the same problem you posted in the maths thread and earlier in this thread. The only requirement we need from the flashlight is that if you shine it orthogonally at a wall, that the irradiation is constant over the illuminated part of the wall, which is approximately true for narrow cones.

Okay, here goes: Instead of projecting the light onto a slanted wall, we create a wall orthogonal to the flashlight's direction and look at the shadow of that slanted wall on the orthogonal wall. Then the amount of light falling onto the slanted wall is equal to the amount of light that would fall onto the shadow of the slanted wall (which is now a quadrilateral taken out of a circle). And since the shadow is orthogonal to the flashlight's direction, it would get constant irradiance, so all you need now is to find the proportion of the light circle that is occupied by the quadrilateral shadow. Essentially all we need to do is intersect a quadrilateral with a circle and find the resulting area. How do we do that?

Easy: Connect all four sides of the quadrilateral to the center of the circle, to split the quadrilateral into four triangles. Four triangles, each with a corner that conveniently lies on the center of the circle. Intersecting such a triangle with that circle will give us a shape that is the sum of at most two circle sectors and at most one triangle. You can calculate the area of circle sectors, you can calculate the area of triangles. So the resulting shape should be no problem for you to calculate the area of:



Note that if the midpoint of the circle is outside of the quadrilateral, you'll have to count some of the triangle areas negatively. I recommend that you use lots of dot and cross products to get all your signs right, and the only trig function you should need is cos^-1. Have fun!

[Mego]

So I've decided to break down and learn PERL. Does anybody have a good tutorial they'd recommend?

[da_nang]

-snip-

I suppose I could apply Green's theorem there?

I also suppose it will be slightly harder. The problem is illustrated by this image:


I need to find the solid angle of the red area (in the xy-plane) subtended at point P, where a cone with an apex angle α (teal) points at direction p (grey), which is a unit vector defined with spherical coordinates as (1,θ,φ) where θ is the longitude and φ is the latitude (both in cyan). The blue conic section is the intersection of the magenta cone with the xy-plane.

So yes, while I know how to project that rectangle to an orthonormal basis of the cone and then project that projection to the plane orthogonal to p, which circle do I use? Assuming I choose the correct one, I can then use Green's theorem to calculate the area.

[MagmaMcFry]

Oh, so you absolutely need the solid angle? Okay, in that case I suggest yet another approach, which still doesn't require any integration, but gets an exact solution anyway. Let's still centrally project the rectangle onto a plane orthogonal to p, so our cone of light intersects the new plane in a circle, not just any ellipse. That's the circle you need to use, by the way. We don't have constant irradiance anymore, but we still have enough symmetry to calculate the exact amount without approximation or integration. Again, let's split the circle and quadrilateral into parts like in my image, and calculate their solid angles individually. Now the solid angle of circle sectors is easily calculated, it's exactly the solid angle of the entire cone of light, multiplied by the fraction of the circle occupied by the circle sector, because symmetry.

To calculate the solid angle of a triangle, we can centrally project the triangle onto an unit sphere around P and calculate the area of the projected triangle. Luckily there are formulas for that. TL;DR just take the sum of the interior spherical angles of the spherical triangle and subtract pi.

There you go, hope that helps.

So I've decided to break down and learn PERL. Does anybody have a good tutorial they'd recommend?

1) It's called Perl. 2) After browsing through some tutorials, I'd say this one is acceptable. Make sure to write lots of tiny programs while reading that tutorial. Also I recommend you skip the database and CGI chapters, unless you need them.

[da_nang]

Let's still centrally project the rectangle onto a plane orthogonal to p, so our cone of light intersects the new plane in a circle, not just any ellipse. That's the circle you need to use, by the way.

Well, see, there an infinite number of planes I could project the rectangle onto, the circle of each being a different size.

If x,y,z are Cartesian coordinates with the positive z-axis being the axis of the cone (I don't care how the other two are oriented as long as their components are orthonormal) and the origin at the apex, then the circle is given by x² + y² = (z tan α)² with z >= 0. So which z am I suppose to use here? Changing z would change the size of the red area when projected onto the plane.

Also, for the spherical trigonometry bit: If the triangle has points P₁, P₂ and P₃, then the normal on the unit sphere where these points are projected are the position vector of each point, yes? As such, the spherical angles would therefore be the angles between these vectors. (I guess I also need to find the edges of the triangle so I can calculate the right angles).

EDIT: Wait, no it wouldn't, but I could calculate the arc lengths and from there use the law of cosines for spherical trigonometry to get them, right?

EDIT2:
Simplest way to find the edges I guess would be to do a Graham scan.

[MagmaMcFry]

Let's still centrally project the rectangle onto a plane orthogonal to p, so our cone of light intersects the new plane in a circle, not just any ellipse. That's the circle you need to use, by the way.

Well, see, there an infinite number of planes I could project the rectangle onto, the circle of each being a different size.

If x,y,z are Cartesian coordinates with the positive z-axis being the axis of the cone (I don't care how the other two are oriented as long as their components are orthonormal) and the origin at the apex, then the circle is given by x² + y² = (z tan α)² with z >= 0. So which z am I suppose to use here? Changing z would change the size of the red area when projected onto the plane.

Pick one of those planes, any one. Each will give the same result. Remember, you're not looking for the area of the red shape, but for its solid angle from P.

Also, for the spherical trigonometry bit: If the triangle has points P₁, P₂ and P₃, then the normal on the unit sphere where these points are projected are the position vector of each point, yes? As such, the spherical angles would therefore be the angles between these vectors. (I guess I also need to find the edges of the triangle so I can calculate the right angles).

EDIT: Wait, no it wouldn't, but I could calculate the arc lengths and from there use the law of cosines for spherical trigonometry to get them, right?

Exactly.

[Rose]

So I have some extracted code here: http://pastebin.com/AfTUwDrk (I think this is everything relevant)

The expected result is

Code: [Select]

Tempout id is 14 and output id is 14(14 can be another number, depending)

Here's what I actually get:

Code: [Select]

Tempout id is 14 and output id is 0
Any idea what I could be doing wrong?

[Bauglir]

Is there any fast way to tell Python that I want comparison operations to treat each lower case letter as being between its upper case and the upper case letter immediately following? IE A<a<B<b...

I could write special behavior at each comparison, but that would bloat the hell out of my code, or I could use the built-in sort function with some extra arguments to get the right behavior whenever I output something, since it's purely a readability thing, but that would be sorting a tree that's ostensibly already sorted, which bugs me irrationally. If there aren't alternatives, or it's profoundly unsafe (I'm gambling on this one being the case), I'll probably go with the latter, but it'd be nice to know.

[Arx]

It may be possible to overload the </> operators, but I'm not certain, and if it is I can't remember the syntax off the top of my head.

Edit: I may or may not have horribly misunderstood your post. If what I said makes no sense contextually, just ignore it.

[Skyrunner]

If it's a sorting function that accepts a predicate, it's totally possible :v

[Levi]

Is there any fast way to tell Python that I want comparison operations to treat each lower case letter as being between its upper case and the upper case letter immediately following? IE A<a<B<b...

I came up with this:

Code: [Select]

list = ['B','D','d','g','C','a','A','c','h']
new = sorted(list, key=lambda e: ord(e.upper()) * 2 + (0 if e.isupper() else 1 ) )
print repr(new)
Edit:  Er, actually I may have misread the question, but you could probably use the code in the lambda to do other comparisons other than sorts.

Maybe this:

Code: [Select]

def letterGT(a,b):
  return ord(a.upper()) * 2 + (0 if a.isupper() else 1 ) > ord(b.upper()) * 2 + (0 if b.isupper() else 1 )

def letterLT(a,b):
  return ord(a.upper()) * 2 + (0 if a.isupper() else 1 ) < ord(b.upper()) * 2 + (0 if b.isupper() else 1 )


[Bauglir]

Ah, yeah. That would amount to writing in a sort of special behavior at each point of comparison. To elaborate, I've got a binary search tree, and want to allow the keys for a given node to be any immutable object. It's just that in the specific case of strings, I don't like the default behavior because of this, just because it means that if you, the user, are reading an output that contains a list of all the nodes, you first have to find where whatever you're looking for ought to be if capitalized, and if it's not there skip to where it would be if it's not.

Although this has actually just given me the idea of forcing it to be all lower case, and implementing some kind of metadata to track where capital letters were prior to the conversion, possibly in the form of an appended number that gets stripped out whenever you print the key for user viewing.

[MagmaMcFry]

I'd really suggest keeping your binary tree sorted with the default comparison, it's much faster that way. You can always reshuffle the keys when you want to display them. Seriously, keeping the objects in viewing order is absolutely not the search tree's responsibility.

[Bauglir]

True, yeah. The thing I just said would slow down insertion more than just re-sorting it whenever it displays, in all likelihood, since I'd have to do a lot of fiddling with the keys to process them.