# Creates fake pulse profile of a given duration (1 sec default)
# that contains a qpo of the given frequency and frequency width
import numpy as np

def generate_pulse_profile(central_frequency,gamma_frequency,output_file,window_duration=None):
    if window_duration is None:
        duration=1
    else:
        duration=window_duration
    npoints=2**12
    ff=np.linspace(0,npoints-1,npoints)
    gamma=gamma_frequency*duration
    f0=central_frequency*duration
    ps=gamma/2/np.pi/((ff-f0)**2+(gamma/2)**2)
    ddtilde=ps**0.5*np.exp((0+1j)*np.random.uniform(0,2*np.pi,npoints))
    dd=np.fft.irfft(ddtilde,npoints)
    dd=dd-np.min(dd)*1.01
    tt=np.linspace(0,duration,npoints)
    np.savetxt(output_file,np.transpose([tt,dd]))
    return

def main():
    import sys
    if len(sys.argv)<4:
        print("""Format:
python fakqpo.py _central_frequency_ _frequency_width_ _Output_File_ [_window_duration_]
""")
    else:
        if len(sys.argv)>4:
            generate_pulse_profile(float(sys.argv[1]),float(sys.argv[2]),sys.argv[3],window_duration=float(sys.argv[4]))
        else:
            generate_pulse_profile(float(sys.argv[1]),float(sys.argv[2]),sys.argv[3])
            
    
# nbins=1000.0
# spec,b=np.histogram(plist, bins=nbins)
# plt.plot(b[1:],spec,'k-')
# plt.plot(chan,cnt*mf*nchan/nbins,'b-')
# plt.show()

if __name__ == "__main__":
    # execute only if run as a script
    main()