//Doubly rotating mirror with lasers simulation
//Assume the primary rotation axis is the vertical = Z+
//The amount of primary rotation is theta
//Since direction of secondary rotation changes sign at every half primary rotation, primary rotation period = 360 deg = 2*pi
//The secondary rotation axis is in the XY plane
//The amount of secondary rotation is phi = mirror angle
//Configured such that as theta increases, phi also increases and increasing phi = reflection moving upwards, increasing theta = reflection moving CCW from above
//If the mirror is double sided, secondary rotation period = 180 deg = pi
//The lasers are in the XY plane arranged uniformly in a circle along theta

//8:68:38->16x, 8:68:42->16x
//8:60:42->5x, 8:60:46->15x, 8:60:38->15x
//94:72:44->~12x at 5.45 Hz nearest to design goal
//92:72:42->~5x with 2x overlaps
//94:72:42->~6x at 8.57 Hz I like this one the most so far
//93:72:42->~4x with overlaps
num_lasers=16;
lasers_theta_offset=0/180*%pi;
fixed_gear_teeth=64;
rotating_gear_teeth=36;
mprintf("Lasers: %g, Fixed teeth: %g, Rotating teeth: %g\n", num_lasers,fixed_gear_teeth, rotating_gear_teeth);
gx_revs=double(lcm([fixed_gear_teeth, rotating_gear_teeth]))/(fixed_gear_teeth);
gx2_revs=double(lcm([fixed_gear_teeth*2, rotating_gear_teeth]))/(2*fixed_gear_teeth);
if(gx_revs~=gx2_revs)
    mprintf("Warning: Return to same position with mirror rotated by half a turn\n");
end
mprintf("Primary revolutions per cycle: %g\n", gx2_revs) //number of primary revolutions to return to identical state (on same side and on either side of mirror)
sd_revs=0.5*(rotating_gear_teeth/fixed_gear_teeth); //number of primary revolutions before secondary mirror is horizontal (either side)
ps_revs=gx2_revs/sd_revs; //number of thetas in which secondary mirror is horizontal prior to return to identical state
mprintf("Secondary half-revolutions per cycle: %g\n", ps_revs)
//lmf=lcm([ps_revs,num_lasers])/num_lasers;
//mprintf("Laser multiplication factor: %g\n", lmf) //not very accurate
ps_thetas=(1:ps_revs)*2*%pi*sd_revs; //theta locations where secondary mirror is horizontal on same (rising or falling) edge (either side of mirror)
ls_thetas=(1:num_lasers)/num_lasers*2*%pi+lasers_theta_offset; //theta locations of lasers

function r=ierem(n,d)
    r=n-d*round(n/d); //IEEERemainder() returns +-distance to nearest multiple
endfunction

function [q,w]=contregion(b)
    //returns vector [a1,a2;c1,c2;...] where each row is start and end index (inclusive) of continuous region of T values in b
    //and boolean w is T if the first and last rows should be combined into a single segment
    a=length(b);
    r=b(1);
    r1=1;
    q=[];
    for i=2:a
        if r==b(i)
            continue;
        end
        if r
            r2=i-1;
            q=[q;r1,r2];
        else
            r1=i;
        end
        r=b(i);
    end
    if r then
        r2=a;
        q=[q;r1,r2];
    end
    w=b(1)&b(a);
endfunction

function plottracelines(rvs,lv,cl)
    [q,w]=contregion(lv);
    if(isempty(q))
        return;
    end
    //if(w)
    //    ts=[tls(q($,1):q($,2)),tls(q(1,1):q(1,2))];
    //    ps=[pls(q($,1):q($,2)),pls(q(1,1):q(1,2))];
    //    q($,:)=[];
    //else
        //ts=tls(q(1,1):q(1,2));
        //ps=pls(q(1,1):q(1,2));
    //end
    rv=rvs(:,q(1,1):q(1,2));
    nq=size(q,1);
    ni=1;
    while ni<=nq
        /*if(fl)
            x=cos(ts).*cos(ps);
            y=sin(ts).*cos(ps);
        else
            x=ts;
            y=ps;
        end*/
        //z=sin(ps);
        param3d(rv(1,:),rv(2,:),rv(3,:));
        //plot(x,y,cl);
        gce().foreground = color(cl);
        ni=ni+1;
        if ni<=nq
            rv=rvs(:,q(ni,1):q(ni,2));
            //ts=tls(q(ni,1):q(ni,2));
            //ps=pls(q(ni,1):q(ni,2));
        end
    end
endfunction

function rv=reflect(sv,mv)
    dv=[mv(1,:)-sv(1);mv(2,:)-sv(2);mv(3,:)-sv(3)];
    dp=sum(dv.*mv,1);
    a1=[mv(1,:).*dp;mv(2,:).*dp;mv(3,:).*dp];
    a2=2.*(dv-a1);
    rv=[sv(1)+a2(1,:);sv(2)+a2(2,:);sv(3)+a2(3,:)];
endfunction

if(%f)
plot(cos(ps_thetas),sin(ps_thetas),'or');
plot(1.1*cos(ls_thetas),1.1*sin(ls_thetas),'o');
isoview on
end
if(%t)
laser_acceptance_angle_h=45/180*%pi; //+-[rad]
laser_acceptance_angle_v=45/180*%pi; //+[rad]
theta_end=gx2_revs*2*%pi;
theta=linspace(0,theta_end,10000);
theta($)=[];
phi=theta*(fixed_gear_teeth/rotating_gear_teeth);
pls=ierem(phi,%pi);
pl1=(pls>=0)&(pls<=laser_acceptance_angle_v);
pl2=(pls<=0)&(pls>=-laser_acceptance_angle_v);
mv=[cos(phi).*cos(theta);cos(phi).*sin(theta);sin(phi)];
mvi=[cos(phi+%pi).*cos(theta);cos(phi+%pi).*sin(theta);sin(phi+%pi)];
clf();
for il=1:num_lasers//ceil(num_lasers/4)
    tls=theta-ls_thetas(il);
    tls=ierem(tls,2*%pi);
    tl1=(abs(tls)<=laser_acceptance_angle_h);
    tl2=(abs(tls)>=(%pi-laser_acceptance_angle_h));
    sv=[cos(ls_thetas(il));sin(ls_thetas(il));0];
    rv=reflect(sv,mv);
    //tsr=ls_thetas(il)+2*tls; //wrong!
    //psr=2*pls;
    lv1=(tl1&pl1); //laser reflected with phi rising
    plottracelines(rv,lv1,'red');
    //rv=reflect(sv,mvi);
    //tsr=ls_thetas(il)+2*(tls-%pi); //wrong!
    //psr=2*-pls;
    lv2=(tl2&pl2); //laser reflected with phi falling
    //plottracelines(rv,lv2,'blue');
end
isoview on
end
