mat_inc_obs=zeros(17,18);
for(incidence=1:1)
    //Calculate wave reflections
    sun_angle=20*incidence;//incidence*5;//degrees up from horizon
    sun_angle=sun_angle/180*%pi;//convert to rad
    wave_height=1.5;//peak to trough
    wave_width=10;//wave_height*3.33 is sharpest wave, anything else becomes flatter
    wave_precision=2000;//use 100 points for the calculation
    wave_sym=wave_precision/2;//50 is a symmetric wave, 2 to 50 has peaks closer to observer, 50 to 98 has peaks closer to horizon
    wave_x=linspace(0,wave_width,wave_precision);//uniform x (initial)
    wave_t=[linspace(0,%pi,wave_sym) linspace(%pi,2*%pi,wave_precision-wave_sym)];//angular time to generate the wave shape (biased)
    wave_x=wave_x+(wave_height/2)*sin(wave_t);
    wave_y=wave_height/2-(wave_height/2)*cos(wave_t);
    //need to find slope at each point (=dy/dx)
    dy=diff(wave_y);
    dx=diff(wave_x);
    surf_angle=atan(dy,dx);//0=flat, +=CCW tilt, -=CW tilt
    ref_angle=sun_angle-2*surf_angle;//reflected angle to observer, 0=horizon, +=up, -=down, 90+=backwards?
    //ref_angle(ref_angle<0)=%nan;//remove negative (=downward) angles
    ref_slope=tan(ref_angle);//slope for collision calculation
    ddx=wave_x-(wave_x(wave_sym)-wave_width);//previous peak we might hit
    ddx(wave_sym:wave_precision)=wave_x(wave_sym:wave_precision)-wave_x(wave_sym);//this peak we might hit
    ddy=wave_y(2:wave_precision)+ref_slope.*ddx(2:wave_precision);//remove one value to use with slope (this is not conservative, ie it will overestimate #transmitted slightly)
    ref_angle(ddy<wave_height)=%nan;//remove any downward or n-reflection angles
    ddx2=linspace(wave_x(wave_sym),wave_width,wave_precision-wave_sym);//now find where incident light may strike based on next wave peak
    ddy2=wave_height-tan(sun_angle)*(wave_x(wave_sym)+wave_width-ddx2);//incident light ray path
    //plot(wave_x,wave_y);
    //plot([0 wave_width/2],[0 wave_width/2*tan(sun_angle)]);
    //plot(wave_x(2:wave_precision),ref_angle/%pi*18);
    if(ddy2(wave_precision-wave_sym)>0)
        [n,m]=min(abs(ddy2-wave_y(wave_sym+1:wave_precision)));
        ref_angle(wave_sym+m:wave_precision-1)=%nan;//light can't strike this area
    end
    for(observ=0:17)
        mat_inc_obs(18-incidence,18-observ)=sum((ref_angle>observ*%pi/36) & (ref_angle<=(observ+1)*%pi/36))
    end
end
//f=scf();
//Matplot(mat_inc_obs);
//f.color_map=hotcolormap(max(max(mat_inc_obs)));
//obs_distr=histc(wave_precision/50,ref_angle/%pi*180,normalization=%f);
//plot([0 wave_width/2],[0 wave_width/2*tan(sun_angle)]);
//plot(wave_x,wave_y);
plot(linspace(0,90,length(ref_angle)),ref_angle/%pi*180);
//plot(wave_x(2:100),ddy)
//plot(ddx2,ddy2);
histplot(wave_precision/100,ref_angle/%pi*180,normalization=%f);
