{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "import numpy as np\n",
    "from IPython import embed\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def GetSASE(CentralEnergy, dE_FWHM, dt_FWHM, samples=0, Axis=True):\n",
    "    h=4.135667662 #in eV*fs\n",
    "    dE=dE_FWHM/2.355 #in eV, converts to sigma\n",
    "    dt=dt_FWHM/2.355 #in fs, converts to sigma\n",
    "    if samples == 0:\n",
    "    \tsamples=int(400.*dt*CentralEnergy/h)\n",
    "    else:\n",
    "    \tif (samples < 400.*dt*CentralEnergy/h):\n",
    "    \t\tprint(\"Number of samples is a little small, proceeding anyway. Got\", samples, \"prefer more than\",400.*dt*CentralEnergy/h)\n",
    "\n",
    "    EnAxis=np.linspace(0.,20.*CentralEnergy,num=samples)\n",
    "    EnInput=np.zeros(samples, dtype=np.complex64)\n",
    "    #for i in range(samples):\n",
    "    EnInput=np.exp(-(EnAxis-CentralEnergy)**2/4./dE**2+2*np.pi*1j*np.random.random(size=samples))\n",
    "    En_FFT=np.fft.fft(EnInput)\n",
    "    TAxis=np.fft.fftfreq(samples,d=(20.*CentralEnergy)/samples)*h\n",
    "    TOutput=np.exp(-TAxis**2/4./dt**2)*En_FFT\n",
    "    EnOutput=np.fft.ifft(TOutput)\n",
    "    if (Axis):\n",
    "    \treturn EnAxis, EnOutput, TAxis, TOutput\n",
    "    else:\n",
    "    \treturn EnOutput, TOutput\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# set the main parameters here:\n",
    "CentralEnergy=80. # in eV\n",
    "bandwidth=0.3 # bandwidth in %\n",
    "dt_FWHM=30.  # FWHM of the temporal duration on average \n",
    "\n",
    "dE_FWHM=CentralEnergy/100 *bandwidth # calculate bandwidth of the spectrum in eV\n",
    "\n",
    "# calculate 3 SASE pulses\n",
    "EnAxis, EnOutput, TAxis, TOutput = GetSASE(CentralEnergy=CentralEnergy, dE_FWHM=dE_FWHM, dt_FWHM=dt_FWHM)\n",
    "EnAxis2, EnOutput2, TAxis2, TOutput2 = GetSASE(CentralEnergy=CentralEnergy, dE_FWHM=dE_FWHM, dt_FWHM=dt_FWHM)\n",
    "EnAxis3, EnOutput3, TAxis3, TOutput3 = GetSASE(CentralEnergy=CentralEnergy, dE_FWHM=dE_FWHM, dt_FWHM=dt_FWHM)\n",
    "\n",
    "\n",
    "# plot them\n",
    "ax1 = plt.subplot(1, 2, 1)\n",
    "plt.plot(EnAxis,np.absolute(EnOutput),EnAxis2,np.absolute(EnOutput2),EnAxis3,np.absolute(EnOutput3) )\n",
    "plt.xlim(CentralEnergy-2.*dE_FWHM,CentralEnergy+2.*dE_FWHM)\n",
    "plt.title('Average pulse duration: %.1f  fs' % dt_FWHM ) \n",
    "ax1.set_xlabel('Photon energy in eV')\n",
    "ax1.set_ylabel('spectral intensity')\n",
    "\n",
    "ax1 =plt.subplot(1, 2, 2)\n",
    "plt.plot(TAxis,np.absolute(TOutput),TAxis2,np.absolute(TOutput2), TAxis, np.max(np.absolute(TOutput))*np.exp(-TAxis**2/2/dt_FWHM**2*2.355**2))\n",
    "plt.plot(TAxis,np.absolute(TOutput),TAxis2,np.absolute(TOutput2), TAxis3,np.absolute(TOutput3))\n",
    "\n",
    "plt.xlim(-2.*dt_FWHM,+2.*dt_FWHM)\n",
    "ax1.set_xlabel('time in fs')\n",
    "ax1.set_ylabel('pulse amplitude')\n",
    "\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}