generate_files.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This script generates all files needed to recreate the data for the manuscript. ~140 GB of storage is needed for the files and ~40 GB of RAM was needed to run certain cells of code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %matplotlib inline\n",
    "import numpy as np\n",
    "import fourier_prop_1D as fp\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import tqdm\n",
    "from scipy.optimize import minimize \n",
    "from scipy.optimize import curve_fit\n",
    "from functools import partialmethod\n",
    "import time as ts\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", message=\"invalid value encountered in arcsin\")\n",
    "\n",
    "# help(fp)\n",
    "\n",
    "#uncomment to disable tqdm output\n",
    "#tqdm.tqdm.__init__ = partialmethod(tqdm.tqdm.__init__, disable=True)\n",
    "\n",
    "#define font sizes used for plotting\n",
    "SMALLER_SIZE = 12\n",
    "SMALL_SIZE = 16\n",
    "MEDIUM_SIZE = 22\n",
    "BIGGER_SIZE = 30\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#the file directory that us used to save the files\n",
    "FILE_DIR='./'\n",
    "\n",
    "#load custom colormap used in plotting\n",
    "cmap = np.load(\"./colormap.npy\")\n",
    "cmap = mpl.colors.ListedColormap(cmap, name='myColorMap', N=cmap.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the laser parameters and simulation window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#center wavelength\n",
    "wvl_0 = 800e-7 #cm\n",
    "\n",
    "#FWHM of the spectrum\n",
    "delta_wvl =100e-7 #cm\n",
    "\n",
    "#1/exp(2) input beam radius of the field intensity\n",
    "w_in = 0.05\n",
    "\n",
    "#length of the spatial simulation window\n",
    "Lx = 8. #cm\n",
    "\n",
    "#number of spatial grid points in simulation\n",
    "Nx = 2**9 + 1\n",
    "\n",
    "#calculate the 1/exp(2) radius of the spectral intensity used in the calculations\n",
    "delta_omega = 2 * np.pi * 2.99792458e10 * delta_wvl / (np.sqrt(2 * np.log(2)) * wvl_0**2) #rad/s\n",
    "\n",
    "#length of the frequency window\n",
    "L_omega = 120.63552 * delta_omega #rad/s\n",
    "\n",
    "#number of frequency grid points in simulation \n",
    "N_omega=  2**13 + 1\n",
    "\n",
    "#focal length of lens\n",
    "focus = 10 #cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define fourier propagation parameters for angle of incidence (AOI): 22.8, 36.3, 52.8 degrees, and the linear approximation\n",
    "aoi228=fp.FourierProp(L=Lx, N=Nx ,N_omega=N_omega, L_omega=L_omega, chirp_pattern='linear_grating_centered', BAR=20, deltax=0., separation=1.97236328, d=1/(1480*10),    aoi=np.deg2rad(22.8), w_in=w_in, focus = focus, delta_wvl=delta_wvl, wvl_0=wvl_0)\n",
    "aoi36=fp.FourierProp(L=Lx, N=Nx, N_omega=N_omega, L_omega=L_omega, chirp_pattern='linear_grating_centered', BAR=20, deltax=0., separation=4.45043945, d=1/(1480*10), aoi=np.deg2rad(36.3), w_in=w_in, focus = focus, delta_wvl=delta_wvl, wvl_0=wvl_0)\n",
    "aoi528=fp.FourierProp(L=Lx, N=Nx, N_omega=N_omega, L_omega=L_omega, chirp_pattern='linear_grating_centered', BAR=20, deltax=0., separation=6.59560547, d=1/(1480*10), aoi=np.deg2rad(52.8), w_in=w_in, focus = focus, delta_wvl=delta_wvl, wvl_0=wvl_0)\n",
    "aoi_linear=fp.FourierProp(L=Lx, N=Nx, N_omega=N_omega, L_omega=L_omega, chirp_pattern='linear', BAR=20.98034668, deltax=0., separation=4.45043945, d=1/(1480*10), aoi=np.deg2rad(36.3), w_in=w_in, focus = focus, delta_wvl=delta_wvl, wvl_0=wvl_0)\n",
    "\n",
    "#store field names and fields. Make sure both arrays are in the same order\n",
    "field_names=['aoi228','aoi36','aoi528','aoi_linear']\n",
    "fields=[aoi228,aoi36,aoi528,aoi_linear]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generates files for the propagation of the fields in (x,z) space at various times. Each file generated is ~35 GB. Also generates files that track the peak intensity througout propagation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n,field in enumerate(fields):\n",
    "    FILE_NAME=field_names[n]\n",
    "\n",
    "    #full length of potentail longtitudinal z-values to propagate through.\n",
    "    #uses the time array of the simulation and converts to z using z=c*t\n",
    "    zs_full=field.timeshift*1e15*0.3/1e4\n",
    "    #use only every 8th term due to memory and storage requirements\n",
    "    zs=zs_full[::8]\n",
    "    \n",
    "    #define arrays for the z_slices of propagation and their local time arrays\n",
    "    z_img=np.zeros((len(zs),field.N_omega,field.N))\n",
    "    z_times=np.zeros((len(zs),field.N_omega))\n",
    "\n",
    "    #construct to field and propagate to focus\n",
    "    field.construct_eField(plot=False)\n",
    "    field.fraunhofer()\n",
    "\n",
    "    #loop through z values and store date in arrays\n",
    "    for i,z in enumerate(tqdm.tqdm(zs)):\n",
    "        field.aspw(zs=np.array(z),plot=False)\n",
    "        z_img[i]=np.abs(field.eFieldxt2[::-1,:])**2\n",
    "        z_times[i]=field.timeshift*1e15+z*1e4/0.3\n",
    "\n",
    "    #save to file\n",
    "    np.save(FILE_DIR+'z_img_'+FILE_NAME, z_img)\n",
    "    np.save(FILE_DIR+'z_times_'+FILE_NAME, z_times)\n",
    "\n",
    "\n",
    "    #extract values from fields\n",
    "    timeshift=field.timeshift #seconds\n",
    "    x2_max=field.x2.max()*1e4 #um\n",
    "\n",
    "    #define arrays\n",
    "\n",
    "    #array to store (x,z) fields that will be plotted\n",
    "    z_plot_t0=np.zeros((len(z_times),Nx))\n",
    "    #also do every 8 time slice due to memory and storage\n",
    "    times=timeshift[::8] #seconds\n",
    "\n",
    "    #define x and z spatial dimensions from simulation\n",
    "    x2=np.linspace(-x2_max,x2_max,Nx) #um\n",
    "    z2=np.linspace(-x2_max,x2_max,len(z_times)) #um\n",
    "    z_max=np.zeros(len(times)) #um\n",
    "    x_max=np.zeros(len(times)) #um\n",
    "    I_max=np.zeros(len(times)) #arb. units\n",
    "\n",
    "    #loop splices together field data that is in (x,t) space into (x,z) space\n",
    "    for i,t in enumerate(tqdm.tqdm(times)):\n",
    "        time_val=np.round(t*1e15,3) #fs\n",
    "\n",
    "        for j in range(len(z_times)):\n",
    "            \n",
    "            if (np.argwhere(np.round(z_times[j],3)==time_val)).size>0:\n",
    "                idx=np.argwhere(np.round(z_times[j],3)==time_val)[0][0]\n",
    "                z_plot_t0[j]=np.array(z_img[j,idx,:])\n",
    "\n",
    "            else:\n",
    "                z_plot_t0[j]=np.zeros((Nx))\n",
    "    \n",
    "        #define where the peak intensity is of the field\n",
    "        max_arg=np.unravel_index((z_plot_t0.T).argmax(), (z_plot_t0.T).shape) \n",
    "        #store location and value of the peak\n",
    "        z_max[i]=(z2[max_arg[1]])\n",
    "        x_max[i]=(x2[max_arg[0]])\n",
    "        I_max[i]=(z_plot_t0.T[max_arg])\n",
    "\n",
    "    #save to file\n",
    "    np.save(FILE_DIR+'z_max_'+FILE_NAME+'.npy',z_max)\n",
    "    np.save(FILE_DIR+'x_max_'+FILE_NAME+'.npy',x_max)\n",
    "    np.save(FILE_DIR+'I_max_'+FILE_NAME+'.npy',I_max)\n",
    "    \n",
    "\n",
    "    #the system I was using was low on memory and so this step was needed to clear the memory before going into another iteration of the loop\n",
    "    #may not be needed on more powerful systsms\n",
    "    del z_img\n",
    "    del z_plot_t0\n",
    "    ts.sleep(15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generates the files containing fits to the peak intensity values of the tracked position to highlight overall trend of data. The fits are qualitative and not meant to represent any specific funcitonal form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#high order polynomial fit\n",
    "def func(x, a, b, c, d, e, f, g, h, i, j):\n",
    "    return a + b*x + c*x**2 + d*x**3 + e*x**4 + f*x**5 + g*x**6 + h*x**7 + i*x**8 + j*x**9\n",
    "\n",
    "#linear fit\n",
    "def func2(x, a, b):\n",
    "    return a + b*x \n",
    "\n",
    "#multiple lorentian fit\n",
    "def func_I(x, a, b, c, d, e, f, g,):\n",
    "    return a + b*np.sqrt(1+x**2/c**2)**-d + e*np.sqrt(1+x**2/f**2)**-g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define the fit for the AOI 22.8 deg case\n",
    "\n",
    "FILE_NAME='aoi228'\n",
    "\n",
    "#load in raw data\n",
    "x_max=np.load(FILE_DIR+'x_max_'+FILE_NAME+'.npy')\n",
    "z_max=np.load(FILE_DIR+'z_max_'+FILE_NAME+'.npy')\n",
    "I_max=np.load(FILE_DIR+'I_max_'+FILE_NAME+'.npy')\n",
    "\n",
    "#define regions to split the fitting to get better results to highlight the overall trend\n",
    "z1=[None]*len(z_max[0:40])\n",
    "x1=[None]*len(x_max[0:40])\n",
    "\n",
    "z2=z_max[40:985]\n",
    "x2=x_max[40:985]\n",
    "\n",
    "z3=[None]*len(z_max[985::])\n",
    "x3=[None]*len(x_max[985::])\n",
    "\n",
    "#define the best fit for the maximum transverse x location\n",
    "popt2, pcov2 = curve_fit(func, z2, x2)\n",
    "\n",
    "\n",
    "#gather data into single array\n",
    "x_fit=np.concatenate((x1,func(z2, *popt2),x3))\n",
    "\n",
    "#define the best fit for the peak intensity\n",
    "popt, pcov = curve_fit(func_I, z_max, I_max, maxfev=60000)\n",
    "\n",
    "#store into array\n",
    "I_fit=func_I(z_max, *popt)\n",
    "\n",
    "#save to file\n",
    "np.save(FILE_DIR+'x_max_'+FILE_NAME+'_fit.npy',x_fit)\n",
    "np.save(FILE_DIR+'I_max_'+FILE_NAME+'_fit.npy',I_fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define the fit for the AOI 36.3 deg case\n",
    "#throws a couple warnings\n",
    "\n",
    "FILE_NAME='aoi36'\n",
    "\n",
    "#load in raw data\n",
    "x_max=np.load(FILE_DIR+'x_max_'+FILE_NAME+'.npy')\n",
    "z_max=np.load(FILE_DIR+'z_max_'+FILE_NAME+'.npy')\n",
    "I_max=np.load(FILE_DIR+'I_max_'+FILE_NAME+'.npy')\n",
    "\n",
    "#define regions to split the fitting to get better results to highlight the overall trend\n",
    "z1=[None]*len(z_max[0:11])\n",
    "x1=[None]*len(x_max[0:11])\n",
    "\n",
    "z2=z_max[11:375]\n",
    "x2=x_max[11:375]\n",
    "\n",
    "z3=z_max[375:650]\n",
    "x3=x_max[375:650]\n",
    "\n",
    "z4=z_max[650:1015]\n",
    "x4=x_max[650:1015]\n",
    "\n",
    "z5=[None]*len(z_max[1015::])\n",
    "x5=[None]*len(x_max[1015::])\n",
    "\n",
    "#define the best fit for the maximum transverse x location\n",
    "popt2, pcov2 = curve_fit(func, z2, x2)\n",
    "popt3, pcov3 = curve_fit(func, z3, x3)\n",
    "popt4, pcov4 = curve_fit(func, z4, x4)\n",
    "\n",
    "\n",
    "#gather data into single array\n",
    "x_fit=np.concatenate((x1,func(z2, *popt2),func(z3, *popt3),func(z4, *popt4),x5))\n",
    "\n",
    "#define the best fit for the peak intensity\n",
    "popt, pcov = curve_fit(func_I, z_max, I_max, maxfev=60000)\n",
    "\n",
    "#store into array\n",
    "I_fit=func_I(z_max, *popt)\n",
    "\n",
    "#save to file\n",
    "np.save(FILE_DIR+'x_max_'+FILE_NAME+'_fit.npy',x_fit)\n",
    "np.save(FILE_DIR+'I_max_'+FILE_NAME+'_fit.npy',I_fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define the fit for the AOI 52.8 deg case\n",
    "\n",
    "FILE_NAME='aoi528'\n",
    "\n",
    "#load in raw data\n",
    "x_max=np.load(FILE_DIR+'x_max_'+FILE_NAME+'.npy')\n",
    "z_max=np.load(FILE_DIR+'z_max_'+FILE_NAME+'.npy')\n",
    "I_max=np.load(FILE_DIR+'I_max_'+FILE_NAME+'.npy')\n",
    "\n",
    "#define regions to split the fitting to get better results to highlight the overall trend\n",
    "z1=z_max[0:510]\n",
    "x1=x_max[0:510]\n",
    "\n",
    "z2=z_max[510:520]\n",
    "x2=x_max[510:520]\n",
    "\n",
    "z3=z_max[520::]\n",
    "x3=x_max[520::]\n",
    "\n",
    "#define the best fit for the maximum transverse x location\n",
    "popt1, pcov1 = curve_fit(func, z1, x1)\n",
    "popt2, pcov2 = curve_fit(func, z2, x2)\n",
    "popt3, pcov3 = curve_fit(func, z3, x3)\n",
    "\n",
    "\n",
    "#gather data into single array\n",
    "x_fit=np.concatenate((func(z1, *popt1),func(z2, *popt2),func(z3, *popt3)))\n",
    "\n",
    "#define the best fit for the peak intensity\n",
    "popt, pcov = curve_fit(func_I, z_max, I_max, maxfev=60000)\n",
    "\n",
    "#store into array\n",
    "I_fit=func_I(z_max, *popt)\n",
    "\n",
    "#save to file\n",
    "np.save(FILE_DIR+'x_max_'+FILE_NAME+'_fit.npy',x_fit)\n",
    "np.save(FILE_DIR+'I_max_'+FILE_NAME+'_fit.npy',I_fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define the fit for the AOI 52.8 deg case\n",
    "\n",
    "FILE_NAME='aoi_linear'\n",
    "\n",
    "#load in raw data\n",
    "x_max=np.load(FILE_DIR+'x_max_'+FILE_NAME+'.npy')\n",
    "z_max=np.load(FILE_DIR+'z_max_'+FILE_NAME+'.npy')\n",
    "I_max=np.load(FILE_DIR+'I_max_'+FILE_NAME+'.npy')\n",
    "\n",
    "#define regions to split the fitting to get better results to highlight the overall trend\n",
    "z1=z_max[0:208]\n",
    "x1=x_max[0:208]\n",
    "\n",
    "z2=z_max[208:210]\n",
    "x2=x_max[208:210]\n",
    "\n",
    "z3=z_max[210:816]\n",
    "x3=x_max[210:816]\n",
    "\n",
    "z4=z_max[816:819]\n",
    "x4=x_max[816:819]\n",
    "\n",
    "z5=z_max[819::]\n",
    "x5=x_max[819::]\n",
    "\n",
    "#define the best fit for the maximum transverse x location\n",
    "popt1, pcov1 = curve_fit(func, z1, x1)\n",
    "popt2, pcov2 = curve_fit(func2, z2, x2)\n",
    "popt3, pcov3 = curve_fit(func, z3, x3)\n",
    "popt4, pcov4 = curve_fit(func2, z4, x4)\n",
    "popt5, pcov5 = curve_fit(func, z5, x5)\n",
    "\n",
    "\n",
    "#gather data into single array\n",
    "x_fit=np.concatenate((func(z1, *popt1),func2(z2, *popt2),func(z3, *popt3),func2(z4, *popt4),func(z5, *popt5)))\n",
    "\n",
    "#define the best fit for the peak intensity\n",
    "popt, pcov = curve_fit(func_I, z_max, I_max, maxfev=60000)\n",
    "\n",
    "#store into array\n",
    "I_fit=func_I(z_max, *popt)\n",
    "\n",
    "#save to file\n",
    "np.save(FILE_DIR+'x_max_'+FILE_NAME+'_fit.npy',x_fit)\n",
    "np.save(FILE_DIR+'I_max_'+FILE_NAME+'_fit.npy',I_fit)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}