CR_Interference.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[]},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["#**Model calculation: Interferometer**\n","\n","> This code supplements Figure 3 of the Chemical Reviews article titled \"Coherence in Chemistry: Foundations and Frontiers\" ([doi.org/10.1021/acs.chemrev.3c00643](https://doi.org/10.1021/acs.chemrev.3c00643))\n","\n","To begin, press the play button beside each cell (proceed sequentially through the cells)."],"metadata":{"id":"H2WPl-rvZSYR"}},{"cell_type":"code","source":["#@title Import libraries, prep plot settings\n","\n","import numpy as np\n","import math\n","import matplotlib.pyplot as plt\n","from scipy import signal\n","from IPython.display import set_matplotlib_formats\n","set_matplotlib_formats('svg')\n","\n","def gensignal(N, trange, Autcrr):\n","    alpha_real = np.random.normal(0, 1/4, int(N))\n","    alpha_im = np.random.normal(0,1/4, int(N))\n","\n","    #initialize empty complex array to store random gaussian samples\n","    alpha = np.zeros(N+1)\n","    alpha = alpha.astype(complex)\n","\n","    for mu in range(1, int(N/2)+1):\n","        real = alpha_real[mu]\n","        im = 1j*alpha_im[mu]\n","        alpha[mu] = real + im\n","        alpha[-mu] = np.conjugate(alpha[mu])\n","\n","    # create alpha for zero freq component\n","    alpha[0] = 1\n","\n","    t= np.linspace(0, trange, N+1)\n","    #construct the spectral representation of the signal A using the random samples alpha_w\n","    fourtransAutcrr = np.fft.fft(Autcrr)\n","    A_w = [alpha[mu]*fourtransAutcrr[mu] for mu in range(-int(N/2), int(N/2)+1)]\n","\n","\n","\n","\n","    A_w = np.fft.ifftshift(A_w)\n","\n","    A_t = np.fft.ifft(A_w)\n","    return np.real(A_t) #discard small imaginary part\n","\n","#calculates the autocorrelation function signal[i]*signal[i+j] averaged over Neave different i's\n","def autocorrelation(signal, j, Neave):\n","    sum = 0\n","    signalj = np.roll(signal, -j)       #creates a copy of the timeseries and shifts it over by j\n","    signal = signal[:Neave]                #remove data from lists that are not included in the ergotic average (Neave is number of timestamps in ergotic average)\n","    signalj= signalj[:Neave]\n","    return np.dot(signal, signalj)/Neave    #returns the ergotic average of signal[i]*signal[i+j] over Neave different i's"],"metadata":{"id":"IksVGVrZzYus","cellView":"form"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#@title Tunable parameters\n","\n","N=2**20\n","dt = 0.01\n","trange = (N+1)*dt\n","t= np.linspace(0, trange, N+1)\n","\n","epsilon = 20\n","tau =   20\n","period = 1\n","omega = 2*math.pi/period"],"metadata":{"id":"-X1cTwAH1cKg"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#@title Calculate and plot\n","\n","\n","###############\n","## Calculate main pulse\n","###############\n","Autcrr = [np.cos(omega*t)*epsilon/tau*np.exp(-t/tau) for t in t]\n","sig = gensignal(N, trange, Autcrr)\n","\n","\n","## Set up delay vectors\n","###################\n","Ndelaysamples = 150 #number of data points in delay plot\n","tdelayrange = 8*period\n","dtdelay = tdelayrange/Ndelaysamples\n","dN = int(dtdelay/dt)\n","t_delay = np.linspace(0, tdelayrange, Ndelaysamples)\n","\n","\n","integrationtime = 2000 * tau\n","Ninteg = int(integrationtime/dt)\n","\n","\n","#Calculate the integrated power on the detector sig[i]sig[i+j] where j is the delay and i sums over the integration time\n","##########\n","detectorsig = np.zeros(Ndelaysamples)\n","\n","for delay in range(0, Ndelaysamples):\n","    sigdelay = np.roll(sig, -delay*dN)\n","    sigdelaydetect = sigdelay[:Ninteg]\n","    sigdetect = sig[:Ninteg]\n","    intensitysum = [(sigdelaydetect[j]+sigdetect[j])**2 for j in range(len(sigdetect))] #instantaneous intensity at time j recorded over the integration time\n","    detectorsig[delay] = np.mean(intensitysum) #average intensity\n","\n","\n","intensitytest = [sig[j]**2 for j in range(len(sig))]\n","val = 2*np.mean(intensitytest)\n","incoherentsum = [val for j in range(Ndelaysamples)]\n","\n","\n","############\n","## Plotting\n","############\n","\n","fig = plt.figure(figsize = (3,6))\n","\n","ax1 = fig.add_subplot(3,1,1)\n","ax1.plot(t, sig, linewidth = 1, color = 'b', label = 'Pulse')\n","ax1.plot(t, np.roll(sig, int(2.4*period/dt)), linewidth = 1, color = 'g', label = 'Delay') #comparison between reference and sample delay\n","plt.xlim(0, 7*period)\n","ax1.axes.xaxis.set_ticks([])\n","ax1.axes.yaxis.set_ticks([])\n","ax1.set_xlabel('Time')\n","ax1.set_ylabel('Signal')\n","ax1.set_title('Pulse Profiles')\n","ax1.legend(prop={'size': 7})\n","\n","ax2 = fig.add_subplot(3,1,2)\n","ax2.plot(t_delay, detectorsig, linewidth = 1, color = 'black', label = 'Interference')\n","ax2.plot(t_delay, incoherentsum, linewidth = 1.5, color = 'blue', linestyle = ':', label = 'Incoherent sum')\n","ax2.set_xlabel('Delay time')\n","ax2.set_ylabel('Average Intensity')\n","ax2.axes.xaxis.set_ticks([])\n","ax2.axes.yaxis.set_ticks([])\n","ax2.legend(prop={'size': 8})\n","ax2.set_title('Detector response')\n","fig.tight_layout()\n","\n","\n","fig2 =  plt.figure(figsize = (1,1))\n","ax3 = fig2.add_subplot(1,1,1)\n","x = np.linspace(0,1, int((N+1)/2)+1)\n","y = np.abs(np.fft.rfft(sig))\n","ax3.plot(x, y, linewidth = 1, color = 'black')\n","plt.xlim(0, 0.05)\n","ax3.axes.xaxis.set_ticks([])\n","ax3.axes.yaxis.set_ticks([])\n","ax3.set_xlabel('Frequency',fontdict={'fontsize': 8})\n","ax3.set_ylabel('Power',fontdict={'fontsize': 8})"],"metadata":{"id":"f932zxF_Wx0b","cellView":"form"},"execution_count":null,"outputs":[]}]}