//Cloud distance circle plot to print out

//input: theta, output width, output height
//units: [rad], [pixel]
//output: line drawn at theta such that a print-out of the equal scaled plot can be held at 45 degrees (one end to top of sky and other end to horizon) and lines will match the theta angle

ths=[95 100 110 120 130 140 150 160 170 175]/180*%pi; //theta angle to draw lines (ellipses)
xm=11; //x extents of output (along 45 degree line)
ym=8.5; //y extents of output (perpendicular to 45 degree line)
a=xm/sqrt(2); //distance param

for(th=ths)

    //ellipse method - not entirely accurate since ellipse center is not 0,0
    //e=th-%pi/4;
    //b=a/sin(e)*sin(%pi-th); //ellipse x-axis
    //d=a/cos(%pi-th)*sin(%pi-th); //ellipse y-axis
    //x=linspace(max(0, b-3/b), min(xm, b),1000);
    //y=(d./b).*sqrt(b.^2-x.^2);
    //x(y>ym+1)=[];
    //y(y>ym+1)=[];
    //plot(x,y,'r');
    
    //projection method - more accurate
    phi=linspace(0, %pi/2, 500);
    t=a./(sin(th-%pi/2)+cos(th-%pi/2).*cos(phi));
    x=cos(th-%pi/2).*cos(phi).*t;
    y=sin(th-%pi/2).*t;
    z=cos(th-%pi/2).*sin(phi).*t;
    b=sign(x).*sqrt(x.^2+(a-y).^2);
    d=z;
    disp([th/%pi*180 b(1)]);
    b(d>ym+1)=[];
    d(d>ym+1)=[];
    plot(b, d);
    
    h = gca(); // get current axes
    h.data_bounds = [0, 0 ; xm, ym];
    h.tight_limits='on';
    h.isoview = 'on';
end
