function [Y21_shot] = ShotNoise(N1, delta)

% Spectrum from Tinto et al.

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

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

for k1 = 1:(N1/2) + 1
    fk = (k1-1)*df; 
    
    Y21_shot(k1) = Y21_shot(k1) * fk * sqrt(18) * (10^(-19)); 
end

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

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