DefaultParams_nonDim

% define two positions to investigate
pointn2_1 = 1.5;
pointp2_1 = 600;

pointn2_2 = 3.5;
pointp2_2 = 600;

parms.n2 = pointn2_1;
parms.p2 = pointp2_1;
N = 1;%5;% 50];
Tend = 12000;
tfinal = Tend;
Output_Times = 0:5:Tend;
totalreps = 1;

mstart = 20*N; % initial conditions
pstart = 333*N;
rstart = 1000*N;

Repeats = 1000;
% Generate average power spec
for i = 1:Repeats
    i
    [ mout,pout,rout ] = GillespieTiming(N,parms,totalreps,mstart, pstart, rstart, tau,Output_Times); 
    pfin = pout(Output_Times>3000);
    fluc = (pfin/N - mean(pfin)/N);
    pf = (N^0.5)*fluc;
    Fs = length(Output_Times)/(tfinal);
    L = length(pf);
    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
    Y = fft(pf,NFFT)/L;
    f = Fs/2*linspace(0,1,NFFT/2+1);
    FourierVec = abs(Y(1:NFFT/2+1));
    f(1) = [];
    FourierVec(1) = [];
    pre = (1/Fs*L)^0.5; 
    FourierVec = pre*FourierVec;
    FourierVec = FourierVec.^2;
    FourierSave(i,:) = FourierVec';
    poutFeedbackLong(i,:) = pfin';
end

AvePower = sum(FourierSave)/Repeats;
AvePower1 = AvePower;
pfin1 = pfin;

%%

parms.n2 = pointn2_2;
parms.p2 = pointp2_2;
N = 1;%5;% 50];
Tend = 12000;
tfinal = Tend;
Output_Times = 0:5:Tend;
totalreps = 1;

mstart = 26*N; % initial conditions
pstart = 68*N;
rstart = 1400*N;

Repeats = 1000;

for i = 1:Repeats
    i
    [ mout,pout,rout ] = GillespieTiming(N,parms,totalreps,mstart, pstart, rstart, tau,Output_Times); 
    pfin = pout(Output_Times>3000);
    fluc = (pfin/N - mean(pfin)/N);
    pf = (N^0.5)*fluc;
    Fs = length(Output_Times)/(tfinal);
    L = length(pf);
    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
    Y = fft(pf,NFFT)/L;
    f = Fs/2*linspace(0,1,NFFT/2+1);
    FourierVec = abs(Y(1:NFFT/2+1));
    f(1) = [];
    FourierVec(1) = [];
    pre = (1/Fs*L)^0.5; 
    FourierVec = pre*FourierVec;
    FourierVec = FourierVec.^2;
    FourierSave(i,:) = FourierVec';
    poutFeedbackLong(i,:) = pfin';
end

AvePower = sum(FourierSave)/Repeats;
AvePower2 = AvePower;
pfin2 = pfin;
%%
% load FIGURE7s2 % to recreate paper fig
figure()
subplot(2,2,1)
plot(2*pi*f,log10(AvePower1))
xlim([0 0.12])
xlabel('Frequency (1/mins)')
ylabel('log_{10}(Power)')
% figure()
a = xlim();
b = ylim();
text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'A'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')

subplot(2,2,2)
plot(Output_Times(Output_Times>3000),pfin1,'color',[0 0.5 0])
xlabel('Time (mins)')
ylabel('Molecule number')
xlim([9000 12000])
a = xlim();
b = ylim();
text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'B'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')


subplot(2,2,3)
plot(2*pi*f,log10(AvePower2))
xlim([0 0.12])
xlabel('Frequency (1/mins)')
ylabel('log_{10}(Power)')
a = xlim();
b = ylim();
text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'C'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')


% figure()
subplot(2,2,4)
plot(Output_Times(Output_Times>3000),pfin2,'color',[0 0.5 0])
xlabel('Time (mins)')
ylabel('Molecule number')
xlim([9000 12000])
a = xlim();
b = ylim();
text(a(1)-0.1*(a(2)-a(1)),b(2)+0.1*(b(2)-b(1)),{'D'},...
    'FontSize',9,'HorizontalAlignment','right', 'VerticalAlignment','bottom','color','k')


