% Looking at how the intensity of the three monitor phosphors vary across time.

% Setting stuff up
boxSize = 100;
Xoffset = 0;
Yoffset = 0;
NHOURS = 6;

disp('Turn on Photometer and hit <return>');
disp('Then focus on white square and hit <return> again');
input('');
PR650init(1,0);

[width, height]=Screen('DisplaySize', 0);
[window,screenRect] = SCREEN(0,'OpenWindow',[190 190 190],[0 0 width height]);
boxRect = [0 0 boxSize boxSize];
boxRect = CenterRect(boxRect,screenRect);
boxRect = boxRect + [Xoffset -1*Yoffset Xoffset -1*Yoffset];
Screen('FillRect',window,[255 255 255],boxRect);
Screen('Flip',window);

input('');
pause(2);
Rspd = []; Gspd = []; Bspd = []; Whitespd = []; timevect = [];
for i = 1:(NHOURS*12)
   Screen('FillRect',window,[255 0 0],boxRect);
   Screen('Flip',window);
   spd = PR650measspd;
   Rspd =[Rspd;spd']
   
   Screen('FillRect',window,[0 255 0],boxRect);
   Screen('Flip',window);
   spd = PR650measspd;
   Gspd =[Gspd;spd'];
   
   Screen('FillRect',window,[0 0 255],boxRect);
   Screen('Flip',window);
   spd = PR650measspd;
   Bspd =[Bspd;spd'];
   
   Screen('FillRect',window,[255 255 255],boxRect);
   Screen('Flip',window);
   spd = PR650measspd;
   Whitespd =[Whitespd;spd'];
   
   Screen('FillRect',window,[190 190 190],boxRect);
   Screen('Flip',window);
   tmp = clock;
   timevect = [timevect; tmp(4) tmp(5)];
   pause(5*60);
end
SCREEN('Close', window);
save TimeCourseOut Rspd Gspd Bspd Whitespd timevect

%%% Analysis
avgRspd = mean(Rspd); 
avgGspd = mean(Gspd); 
avgBspd = mean(Bspd);
avgWspd = mean(Whitespd); 

figure; hold on;
plot(avgRspd,'r-');
plot(avgGspd,'g-');
plot(avgBspd,'b-');
plot(avgWspd,'k-');

tmp = diff(timevect(:,1));
daychanges = find(tmp == -23); % Day changes
for i = 1:length(daychanges)
   timevect([daychanges(i)+1:end],1) =...
   timevect([daychanges(i)+1:end],1) + 24;
end

figure; hold on;
weights = [];
for i = 1:size(Rspd,1)
   weights = [weights regress(avgRspd',Rspd(i,:)')];
end
plot(timevect(:,1)+timevect(:,2)*(1/60),weights,'r-');
weights = [];
for i = 1:size(Gspd,1)
   weights = [weights regress(avgGspd',Gspd(i,:)')];
end
plot(timevect(:,1)+timevect(:,2)*(1/60),weights,'g-');
weights = [];
for i = 1:size(Bspd,1)
   weights = [weights regress(avgBspd',Bspd(i,:)')];
end
plot(timevect(:,1)+timevect(:,2)*(1/60),weights,'b-');
%weights = [];
%for i = 1:size(Whitespd,1)
%   weights = [weights regress(avgWspd',Whitespd(i,:)')];
%end
%plot(timevect(:,1)+timevect(:,2)*(1/60),weights,'k-');
ylabel('Weight to get to average');
xlabel('Time (hours)');


%%% Test of linearity for white
pred = avgRspd+avgGspd+avgBspd;
figure; axes; hold on;
plot(pred,'b-');
plot(avgWspd,'k-');
[weights] = regress(avgWspd',[avgRspd' avgGspd' avgBspd']);
gunnames = {'red','green','blue'};
for i = 1:3
    if weights(i) < 1
        percentchange = 100*(1-weights(i));
        disp([char(gunnames(i)),' gun is ',num2str(percentchange,'%0.2g'),'% too dim in white']);
    else
        percentchange = 100*(weights(i)-1);
        disp([char(gunnames(i)),' gun is ',num2str(percentchange,'%0.2g'),'% too bright in white']);
    end
end