//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
//If the mirror is triple sided, secondary rotation period = 120 deg = pi*2/3
//The lasers are in the XY plane arranged uniformly in a circle along theta

//lasers,mirrors,fixed,rotating (following at +-45deg acceptance)
//94,1,72,42 = 3x multiplication, 12 acc wind, 7 primary revs per cycle, continuous, no overlaps
//94,2,72,42 = 6x multiplication, 24 acc wind, 7 primary revs per cycle, continuous, no overlaps
//98,2,72,42 = 6x multiplication, 24 acc wind, 7 primary revs per cycle, continuous, no overlaps
//94,3,72,42 = 9x multiplication, 36 acc wind, 7 primary revs per cycle, continuous, no overlaps
//94,2,64,40 = 4x multiplication, 16 acc wind, 5 primary revs per cycle, continuous, no overlaps
//94,3,64,40 = 6x multiplication, 24 acc wind, 5 primary revs per cycle, continuous, no overlaps -this could be a reasonable choice if higher resolution is desired along with no mirror degeneracy
//98,3,64,40 = 6x multiplication, 24 acc wind, 5 primary revs per cycle, continuous, no overlaps
//106,3,64,40 = 6x multiplication, 24 acc wind, 5 primary revs per cycle, continuous, no overlaps
//118,3,64,40 = 6x multiplication, 24 acc wind, 5 primary revs per cycle, continuous, no overlaps
//94,3,80,50 = 6x multiplication, 24 acc wind, 5 primary revs per cycle, continuous, no overlaps
//*** 94,3,64,36 = 4x multiplication, 16 acc wind, 3 primary revs per cycle with 3x mirror degeneracy, continuous, no overlaps *presently preferred configuration due to fast update rate and similar xy resolution assuming 1 MHz update clock over 96 bits
//96,3,80,36 = 5x multiplication, 20 acc wind, 3 primary revs per cycle with 3x mirror degeneracy, continuous, no overlaps
//98,3,64,36 = 4x multiplication, 16 acc wind, 3 primary revs per cycle with 3x mirror degeneracy, continuous, no overlap
//92,4,64,40 at +-22.5deg acceptance = 4x multiplication, 32 acc wind, 5 primary revs per cycle, continuous, no overlaps
//92,4,96,36 at +-22.5deg acceptance = 4x multiplication, 32 acc wind, 3 primary revs per cycle, cno
//93,4,72,40 at +-30deg acceptance = 6x multiplication, 36 acc wind, 5 primary revs per cycle, continuous, no overlaps
//93,4,60,40 at +-30deg acceptance = 1x multiplication, 6 acc wind, 1 primary revs per cycle, continuous, no overlaps
//93,4,66,44 at +-30deg acceptance = 1x multiplication, 6 acc wind, 1 primary revs per cycle, continuous, no overlaps
//93,4,60,50 at +-30deg acceptance = 4x multiplication, 24 acc wind, 5 primary revs per cycle, continuous, no overlaps
//96,3,70,40 at +-60deg acceptance = 7x multiplication, 21 acc wind, 4 primary revs per cycle, cno

//Rule for CNO: gcd(2*L,W)=C where L=num_lasers, W=acc_wind, C=360deg/(acceptance)

num_lasers=92;
lasers_theta_offset=0/180*%pi;
num_mirrors=4;
fixed_gear_teeth=96;
rotating_gear_teeth=36;
laser_acceptance_angle_h=22.5/180*%pi; //+-[rad]
laser_acceptance_angle_v=45/180*%pi; //+[rad]

mprintf("NL %g, NM %g, FT %g, RT %g\n", num_lasers,num_mirrors,fixed_gear_teeth, rotating_gear_teeth);
gx_revs=double(lcm([fixed_gear_teeth, rotating_gear_teeth]))/(fixed_gear_teeth);
gxn_revs=double(lcm([fixed_gear_teeth*num_mirrors, rotating_gear_teeth]))/(num_mirrors*fixed_gear_teeth);
if(gx_revs~=gxn_revs)
    mprintf("Warning: Return to same position with mirror rotated to another face\n");
end
mprintf("Primary revolutions per cycle: %g\n", gxn_revs) //number of primary revolutions to return to identical state (on any side of mirror)
sd_revs=(rotating_gear_teeth/fixed_gear_teeth)/num_mirrors; //number of primary revolutions before secondary mirror is horizontal and rising (any side)
ps_revs=gxn_revs/sd_revs; //number of thetas in which secondary mirror is horizontal prior to return to identical state
mprintf("Acceptance windows per cycle: %g\n", ps_revs)

//return

ps_thetas_rising=(0:(ps_revs-1))*2*%pi*sd_revs; //[rad] unwrapped theta locations where secondary mirror is horizontal on rising edge (any side of mirror)
if(modulo(num_mirrors,2)<0.1)
    ps_thetas_falling=ps_thetas_rising-%pi;
else
    ps_thetas_falling=ps_thetas_rising-(0.5*2*%pi*sd_revs)-%pi; //[rad] unwrapped theta locations where secondary mirror is horizontal on falling edge (any side of mirror)
end

ls_thetas=(0:(num_lasers-1))/num_lasers*2*%pi+lasers_theta_offset; //theta locations of lasers

//plot acceptance windows
//find acceptance window center boundaries in theta and phi
ts=ps_thetas_rising(1);
psl=fixed_gear_teeth/rotating_gear_teeth; //rise of phi vs theta for rising edge, also slope of laser projection line dphi/dtheta
dphi=laser_acceptance_angle_v;
dtheta=dphi/psl; //[rad] revolution of theta for one phi acceptance window
if(abs(dtheta+2*laser_acceptance_angle_h)>=(2*%pi))
    error("Acceptance window is more than 2*%pi, cannot continue");
end

function plot_wrap(lproj,col)
    for(i=1:size(lproj,2))
        ts=lproj(1,i);
        td=lproj(3,i)-ts;
        if(td<-1e-6)
            error("plot_wrap does not support negative dtheta");
        end
        if(td<1e-6)
            continue;
        end
        ts=modulo(ts,2*%pi);
        if(ts<0)
            ts=ts+2*%pi;
        end
        te=ts+td;
        ps=lproj(2,i);
        pe=lproj(4,i);
        pd=pe-ps;
        wp=%f;
        while(te>2*%pi)
            pc=(2*%pi-ts)/td;
            plot([ts,2*%pi],[ps,ps+pd*pc],'color',col);
            ps=ps+pd*pc;
            te=te-2*%pi;
            ts=0;
            wp=%t;
        end
        if(wp)
            plot([ts,te],[ps,pe],'color',col);
        else
            plot([ts,te],[ps,pe],'color',col);
        end
    end
endfunction

//awc=[ps_thetas_falling(1)-dtheta,ps_thetas_falling(1);dphi,0]; //window center for falling edge
clf();

cmap=jetcolormap(length(ps_thetas_rising));

for wnum=1:length(ps_thetas_rising)
awc=[ps_thetas_rising(wnum),ps_thetas_rising(wnum)+dtheta,0,dphi]; //window centers for rising edge
aws=[awc(1)-laser_acceptance_angle_h,awc(1)+laser_acceptance_angle_h,awc(2)-laser_acceptance_angle_h,awc(2)+laser_acceptance_angle_h];//where first and last accepted laser should be
dlt=2*%pi/num_lasers;//laser dtheta
ln1=ceil((aws(1)-lasers_theta_offset)/dlt)*dlt;//[rad] unwrapped where first accepted laser actually is
luw=ls_thetas+ln1;//[rad] unwrapped lasers from start of acceptance window
//luw(luw<aws(1))=[];
luw(luw>aws(4))=[]; //keep only lasers inside acceptance window
lproj=zeros(4,length(luw)); //laser projected line endpoints (theta1;phi1;theta2;phi2)
for(i=1:length(luw))
    if(luw(i)>aws(2)) //laser starts late
        lpcs=(luw(i)-aws(2))/dtheta;
    else
        lpcs=0;
    end
    if(luw(i)<aws(3)) //laser ends early
        lpce=1-(aws(3)-luw(i))/dtheta;
    else
        lpce=1;
    end
    awl=awc(1)+dtheta*lpcs;//theta center at start
    lproj(1,i)=awl-(luw(i)-awl);//reflected projected theta
    awp=0+dphi*lpcs;//phi center at start
    lproj(2,i)=awp-(0-awp);//reflected projected phi
    awl=awc(1)+dtheta*lpce;//theta center at end
    lproj(3,i)=awl-(luw(i)-awl);//reflected projected theta
    awp=0+dphi*lpce;//phi center at end
    lproj(4,i)=awp-(0-awp);//reflected projected phi
end
if(%t)
    //plot(aws([1,3,4,2,1]),awc([3,4,4,3,3]),'k'); //acceptance window outline
    //plot(lproj([1,3],:),lproj([2,4],:),'r'); //projected laser locations
    plot_wrap(lproj,cmap(wnum,:));
end
end
title(sprintf("%g acceptance windows",wnum));
if(%f)
    clf();
    plot(cos(ps_thetas_rising),sin(ps_thetas_rising),'or');
    plot(cos(ps_thetas_falling),sin(ps_thetas_falling),'+g');
    plot(1.1*cos(ls_thetas),1.1*sin(ls_thetas),'o');
    isoview on
end

return

if(%f)
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);
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));
    tsr=ls_thetas(il)+2*tls;
    psr=2*pls;
    lv1=(tl1&pl1); //laser reflected with phi rising
    plottracelines(tsr,psr,lv1,'red',%t);
    tsr=ls_thetas(il)+2*(tls-%pi);
    psr=2*-pls;
    lv2=(tl2&pl2); //laser reflected with phi falling
    plottracelines(tsr,psr,lv2,'blue',%t);
end
isoview on
end
