function [ractive rinactive output]= RedLightSimulation()


n=20;
m=20;
rdurv=linspace(0,10000,n);
rintv=linspace(0,10,m);
[X,Y]=meshgrid(rdurv,rintv);

ractive=zeros(n,m);
rinactive=zeros(n,m);
output=zeros(n,m);

global p

for j=1:n
    dur=rdurv(j);
    for k=1:m
        int=rintv(k);
        %display(['----------------------------------------------------']);
        %display(['--------- Emission Duration ' num2str(dur) ' seconds ---------']);
        %display(['----------------------------------------------------']);


        % 1 min 60
        % 1 h 3600
        % 1 d 86400 

        p=struct();
        %red light
        p.rstart=0;
        p.rdur=dur;
        p.rint=int;
        [t,y]=RedLightSensor();

        %% RedLight
        mRNAthresh=0.5;
        rstart=find(y(:,8)>mRNAthresh,1);
        rtstart=t(rstart);
        rtend=t(find(y(rstart+10:end,8)<mRNAthresh,1));
        if size(rtstart-p.rstart,1)>0
            ractive(j,k)=rtstart-p.rstart;
        else
            ractive(j,k)=-20;
        end
        if size(rtend-p.rstart-p.rdur,1)>0
            rinactive(j,k)=rtend-p.rstart-p.rdur;
        else
            rinactive(j,k)=-20;
        end
        %display(['REDLIGHT Time for activation: ' num2str(rtstart-p.rstart) ' s - Time for deactivation: ' num2str(rtend-p.rstart-p.rdur) ' s '])



        %% Output
        output(j,k)=y(end,10);
        %display(['Final output of the dye is ' num2str(output) ' nM'])
        
    end
end

figure(1)
clf
surf(X',Y',ractive)
title('Activation Time Red Light')
xlabel('red exposure time (s)')
ylabel('red intensity')
zlabel('activation time (s)')

figure(2)
clf
surf(X',Y',rinactive)
title('Dectivation Time Red Light')
xlabel('red exposure time (s)')
ylabel('red intensity')
zlabel('deactivation time (s)')

figure(5)
clf
surf(X',Y',output)
title('Output of dye ')
xlabel('red exposure time (s)')
ylabel('red intensity')
zlabel('dye output (nM)')

end