function [lnoise] = LaserNoise(N1, delta)

% Spectrum from Tinto et al and Bender et al (Pre Phase Report)

lnoise = randn(N1, 1); 
lnoise = fft (lnoise);

df = (1/(N1*delta)); 

for k1 = 2:(N1/2) + 1
    fk = (k1-1)*df; 
    
    tf = ( (10^(-28)* ((fk)^(-2/3))) + (6.3* (10^(-37)) * ( (fk)^(-3.4)))); 
    tf = sqrt(tf); 
    
    lnoise(k1) = tf * lnoise(k1);
end

lnoise(1) = lnoise(1)*( sqrt( (10^(-28)) + (6.3*(10^(-37))))) ; 

for k1 = (N1/2)+2: N1
    lnoise(k1) = conj(lnoise( N1 - k1 +2)); 
end

%imaginary part in the order of 10^-34, 
%real part 10^-19. imaginary part : computational artifact
lnoise = real(ifft(lnoise));