## Calculate Relative intensity noise using Rate equations

Hello Everybody,

I'm trying to plot the frequency spectrum of RIN using Matlab.

I've the following code

function p=rateblock
tau_s = 3e-9;
N0 = 1e24;
A =1e-12;
P0 = 1/(A*tau_s);

TSPAN = 0:1e-2:10.23;;
Y0 =[0 0];

[T,Y] = ODE45(@rate_equation,TSPAN,Y0);
subplot(2,1,1)
plot(T*tau_s ,Y(:,1)*N0)
title('carriers density in high laser level')
subplot(2,1,2)
plot(T*tau_s ,Y(:,2)*P0)
title('photons density in activer region') %

dt=1e-2;
s_avg=sum(Y(:,2))/1024;%average photon numbers
ds=Y(:,2)-s_avg;%photon fluctuation with respect to average
TT=1e-6; %total time of integration
N=1024;%fft length
F=[0:N-1]';
RIN=((dt^2)/((s_avg)^2*TT)).*(abs(fft(ds))).^2;
RIN=10*log10(RIN);
plot(F,RIN)
title('RIN')

And for the Rate equation i used this function
function dy = rate_equation(t,y)
dy = zeros(2,1);

tau_s = 3e-9; % carriers lifetime
tau_p = 1e-12; % photons lifetime
A = 1e-12; % linear gain costant
N0 = 1e24; % trasparency carries density
V = 3.75e-14; % modal volume
gamma = 1e-5; % gain compression factor
q = 1.6e-19; % electron charge

I0 = N0*q*V/tau_s; % trasparency current
tau_norm = tau_s/tau_p;
eta = A*tau_p*N0; % efficiency

I = 2.5*I0; % pumping current

dy(1)= I/I0 -y(2)*(y(1) - 1) -y(1);
dy(2) = tau_norm*(y(2)*(eta*(y(1) - 1) -1) + gamma*eta*y(1))

Although I'm facing problems with caluclating the RIN and plot it with frequency ...

Can anybody help me with this issue?