# Pulsar specifications: B0329+54 @ 400MHz
flux= 1.5; # Jansky
period=0.714520; # seconds
width=0.0066; # seconds

# Simulation specs
runs=7;
runtime=period*runs;
samplerate = 200000;

samples=runtime*samplerate;
time = (0:1/samplerate:(samples-1)/samplerate);

# Telescope specifications
Tsys=50; # *System noise figure, in Kelvin
antenna = 0.09; # Kelvin per Janksy

# Simulate the signal as received
phase=rand*period; # Pick a random phase in the pulsar sequence to start

# Recatangular pulse
#pulsar=sqrt(flux*antenna*period/width)*randn(1,samples) .* (sign(width-mod(time+phase,period))+1)/2;

# Triangular pulse with 50% width at W50 value
pulsar=sqrt(flux*antenna*period/width)*randn(1,samples) .* sqrt(max(0,abs(mod(time+phase,period)-period/2)-period/2+width)/width);

noise=sqrt(Tsys)*randn(1,samples); # System noise
input=pulsar+noise;

# Simulate the receiver

# Square law detector
square=input.^2; 

# Integrator: average a block of samples
block = floor(samplerate*width); #Number of samples to integrate over
gain =sqrt(block);
printf("Integrator gain: %d\n",gain)

n=8; # Number of overlapping integrations

for i=0:floor(n*samples/block)-n;
avg(i+1)=sum(square(floor(1+i*block/n):floor((i+n)*block/n)))/block;
endfor;

time2=0:block/samplerate/n:(samples-block)/samplerate;

plot(time2,avg,";Tsys=50K, BW=100kHz, 6.6ms integration, B0329+54 (simulated);")
