In [2]:

    

from __future__ import division
import numpy as np
import matplotlib.pylab as plt
import time

fmin = 311.25*10**6 #Hz
fban = 16*10**6 #Hz
freq = fmin + 0 * fban
fsample = 2*fban
Period = 1.6*10**(-3) #seconds
PulseWidth = Period/200


fref = fmin#Hz

FileName = "data/test2_caustic_"

# This generates the signal time series. It's some gaussian envelope * random values between +-0.5
def Signal(width=PulseWidth*fsample, length=int(Period*fsample), Noise=0):
    s = np.random.rand(length)-0.5
    t = np.array(range(length))
    t0 = length/2
    envelop = np.exp( -((t-t0)/width)**2 )
    return envelop*s

def PhaseArray(SignalBand, GeoPath, DisPath, Freq):
    LG = len(GeoPath)
    LD = len(DisPath)
    if LG > LD:
        DisPath = np.pad(DisPath,int((LG-LD)/2)+1,'edge')
        DisPath = DisPath[:LG]
    w = Freq + np.array(range(SignalBand))*fban/SignalBand # frequency array, of size fourier transform of signal
    PA = np.outer(w,GeoPath) - np.outer(((fref**2)/w),DisPath)  # the "phase array". it consist of a vertical array of frequencies, each of which has some paths associated.
    return PA

def PhaseFactor(PV):
    phase = np.exp(2*np.pi*(0+1j)*PV )
    return phase

def PathInt(PA): # semi-analytical
    PathInt = []
    fRange, pathRange = np.shape(PA)
    for i in range(fRange):
        Onef = 0
        IntRange, weight = FindIntRange(PA[i,:])
        for j in IntRange:
            Onef += weight[j] * PhaseFactor(PA[i,j])
        PathInt += [Onef]
    return np.array(PathInt)

def FindIntRange(TotPath):
    width = 0.289
    dPath = np.gradient(TotPath)
    window = np.exp(-dPath**2/(2*(width)**2))
    IntRange = np.where(window > 3e-3)[0]
    return IntRange, window

def Scan(begin, end, freq):
    scan = range(begin,end)
    res = np.linspace(-1/2,1/2,10,endpoint='true')[:-1]
    lensed = []
    spec = []
    for i in scan:
        dp = dispath[i-(dens*100):i+(dens*100+1)]
        for j in res:
            gp = GeoPath(j)
            PA = PhaseArray(l,gp,dp,freq)
            PI = PathInt(PA)
            s1 = np.fft.irfft(sf*PI)
            # save intensity as well as spectrum
            lensed += [(s1**2).sum()]
            spec += [sf*PI]
        print freq/10**6, i, begin, end, time.clock()
    return np.array(lensed)/norm, np.array(spec)

def DisPath():
    x = np.random.rand(10)*2-1
    y = np.random.rand(10)*2-1
    z = (1+0j)*x + (0+1j)*y
    q = np.array([0]*991)
    z = np.concatenate((z,q))
    z = np.fft.irfft(z)
    amp = 4*10**(-6)
    z = z*amp/np.max(np.abs(z))
    return z


    

In [3]:

    

%cd /home/flin/scr/binary-lens


    

    




/mnt/scratch-lustre/flin/binary-lens

In [4]:

    

%matplotlib notebook


    

In [5]:

    

dispath = np.load(FileName+'Dis.npy')
geopath = np.load(FileName+'Geo.npy')


    

In [6]:

    

plt.figure()
plt.plot(dispath)
plt.plot(geopath)


    

    















    












    
Out[6]:





[<matplotlib.lines.Line2D at 0x7f7ffa1a3c90>]

In [110]:

    

plt.figure()
plt.plot(np.gradient(np.gradient(geopath)))
plt.plot(np.gradient(np.gradient(dispath)))
print len(dispath)


    

    















    












    




2400

In [184]:

    

def DisPath():
    x = np.random.rand(10)*2-1
    y = np.random.rand(10)*2-1
    z = (1+0j)*x + (0+1j)*y
    q = np.array([0]*991)
    z = np.concatenate((z,q))
    z = np.fft.irfft(z)
    amp = 4*10**(-6)
    z = z*amp/np.max(np.abs(z))
    return z

dis = DisPath()

dens = 4
temp = np.fft.rfft(dis)
temp = np.concatenate((temp,np.zeros( (dens-1)*1000)))
temp = np.fft.irfft(temp)
temp *= dens
dis = temp

plt.figure()
plt.plot(np.arange(len(dis)),np.ones(len(dis))*np.amax(np.gradient(np.gradient(geopath))))
plt.plot(np.gradient(np.gradient(dis)))


    

    















    












    
Out[184]:





[<matplotlib.lines.Line2D at 0x7f7de651b7d0>]

In [185]:

    

loc = 2331
loc2= 1200 -400
plt.figure()
plt.plot(dis[loc-1200:loc+1200])
plt.plot(np.arange(loc2,loc2+len(geopath)),geopath)
plt.plot(np.arange(loc2,loc2+len(geopath)),geopath-dis[loc-1200:loc+1200][loc2:loc2+len(geopath)])
print len(dis[loc-1200:loc+1200])


    

    















    












    




2400

In [186]:

    

np.save('data/causticDis15.npy',dis[loc-1200:loc+1200])


    

In [ ]: