\documentclass[pra,onecolumn,reprint]{revtex4}
%使用 revTeX 的 onecolume 版本
\usepackage{amsmath}    
\usepackage{graphicx}  
\usepackage{natbib}  

%\usepackage{epstopdf}
%\pagestyle{empty} %去除頁碼


\begin{document}



\title{Advanced Laser Technologies  Homework \#2}
\author{Chih-Han Lin~~~student ID: r99245002}
\affiliation{Graduate Institute of Applied Physics, National Taiwan University}


\maketitle


\section{Demonstration of mode-locking pulses}
Assume that complex amplitude has following expression,

\begin{align*}
A(t)=\sum_a A_q \exp\left(\frac{iq2\pi t}{T_f}\right)
\end{align*}

\noindent and pulse intensity $I(t)=|A(t)|^2$. We choose modes $M=11$ ($q$ sums from 1 to 11) with different complex coefficient 
$A_q$ as follows,

\begin{itemize}
\item [(a)] Equal magnitude and equal phases, i.e. $A_q=1$.
\item [(b)] Magnitude that obey the Gaussian spectral profile and  equal phase, i.e. $A_q=\exp\left[-\frac{1}{2}\left(\frac{q}{5}\right)^2\right]$.
\item [(c)] Equal magnitude and random phases (uniform distributing between 0 and $2\pi$), i.e. $A_q=\exp(rand()*2\pi)$.
\item [(d)] Equal magnitude and phase that obey the Gaussian distribution, i.e. $A_q=\exp \left[-i\cdot 2\pi\frac{1}{2}\left(\frac{q}{5}\right)^2\right]$
\end{itemize}

\section{Simulation result}
\subsection{Equal magnitude with equal phases}

FIG. 1. shows additive intensity function with equal relative phase in time domain and frequency domain respectively. Square Mask in  frequency domain would generate ripples between pulses in time domain (oscillation with decay property of sinc function).   

\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]{equal-amp-equal-phase.pdf}
\caption{Equal magnitude and equal phases}
\end{figure}

\subsection{Gauassian distributed magnitude with equal phases}

First, we demonstrate Gaussian distributed magnitude with equal phase, i.e. $A_q=\exp\left[-\frac{1}{2}\left(\frac{q-6}{5}\right)^2\right]$, and the simulation result is shown in FIG. 2. Since the square mask effect is almost neglected comparing to standard deviation of Gaussian distributed spectral profile, the additive ultrashort pulse in time domain would have less ripple than previous case (modes with equal amplitude).

\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]{guassian-amp-equal-phase.pdf}
\caption{Gauassian distributed magnitude with equal phases}
\end{figure}

We than increase the mode number ($M$=101) with symmetrical spectrum and observe the dependence between pulse duration and Gaussian spectral width by varying the complex mode's amplitude with:

\begin{align*}
A_q=\sum_a A_q \exp\left(-4\frac{q-50}{50k}\right),~~~\mbox{k}=1,2,\cdots ,5.
\end{align*}

The result shown in FIG. 3. says that pulse duration becomes shorter as $k$ increases. If the standard deviation $\sigma$  of the Gaussian spectrum is large compared with total mode spacing, the shortest pulse duration in time domain is limited by total mode spacing in frequency domain. We may slowly change the $\sigma$ via computer simulation (seeing FIG. 4.) to specify different operation regime of spectrum control.

\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]{guassian-amp-equal-phase-multi.pdf}
\caption{pulse duration decrease with $\sigma$}
\end{figure}      

\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]{pulse-duration-sigma.pdf}
\caption{pulse duration v.s. $\sigma$}
\end{figure}   

\subsection{Equal magnitude with random phases}

If mode relative phase is a random number between $0$ and $2\pi$, the additive pulse would have lower AC component and generate randomly periodic signal which not specified as any usable pulses. The simulation result is shown in FIG. 5.

\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]{equal-amp-rand-phase.pdf}
\caption{equal magnitude with random phases}
\end{figure} 

\subsection{Equal magnitude with Gaussian distributed phases}

We may simulate Gaussian shape relative phase with parameters taking values  $A_q=\exp \left[-i\cdot \theta_c\frac{1}{2}\left(\frac{q}{5}\right)^2\right]$ with $\theta_c=0.5\pi, \pi$ and $2\pi$. We observe the pulse duration broaden as neighbor mode relative phase getting larger (seeing FIG. 6.).

  
\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{equal-amp-guassian-phase.pdf}
\caption{equal magnitude with Gaussian distributed phases}
\end{figure} 

\appendix
\section{Matlab code}

We use a for loop to directly sum the additive pulse in time domain and adjust parameter {\verb amp } and {\verb phase } to simulate different cases (Gaussian distribution or constant).

\begin{verbatim}
T_f=1;
samp_time=1E-4;
t=0:samp_time:3;
E_total=zeros(size(t));
for mode=1:11,
    amp=1;
    %amp=exp(-0.5*(mode/5)^2);
    phase=-0.5*pi*(mode/5)^2;
    %phase=rand()*2*pi;
    %phase=0;
    A_q=amp*exp(j*phase);
    E_total=E_total+A_q*exp(j*mode*2*pi*t/T_f);
end
\end{verbatim}

To reconstruct the spectrum information, we use embedded {\verb fft } function in Matlab.
 

\begin{verbatim}
samp=size(t,2);
Spectrum=fft(E_total,samp);
f=1/samp_time*(1:samp)/samp;
I_f=Spectrum.*conj(Spectrum)/samp;
I_fphase=atan(imag(Spectrum)./real(Spectrum));
\end{verbatim}

Finally we plot intensity function with phase both in time and frequency domain via following codes:

\begin{verbatim}
subplot(221)
plot(t,abs(E_total).^2);
xlabel('time (arbitary unit)');
ylabel('Intensity (arbitary unit)')
subplot(223)
plot(t,atan(imag(E_total)./real(E_total)),'r')
xlabel('time (arbitary unit)');
ylabel('phase (radian)')
subplot(222)
plot(f(1:cutoff_freq),I_f(1:cutoff_freq));
xlabel('Frequency (arbitary unit)');
ylabel('Intensity (arbitary unit)')
subplot(224)
plot(f(1:cutoff_freq),I_fphase(1:cutoff_freq),'r');
xlabel('Frequency (arbitary unit)');
ylabel('phase (radian)')
\end{verbatim}

To generate data like FIG. 3., we have to add some analysis to automatically recognize the FWHM of the short pulse, the modified code shows as follows,

\begin{verbatim}
\end{verbatim}


\begin{verbatim}

T_f=1;
samp_time=1E-4;
t=0:samp_time:3;
num=5;
E_total=zeros(num,size(t,2));
for k=1:num,
    for mode=1:101,
    %amp=1;
    %amp=exp(-4*((mode-50)/(50*k))^2);
    phase=-2*pi*(mode/5)^2;
    %phase=rand()*2*pi;
    phase=0;
    A_q=amp*exp(j*phase);
    E_total(k,:)=E_total(k,:)+A_q*exp(j*mode*2*pi*t/T_f);
    end
a0=0;
It=abs(E_total(k,:)).^2;
for i=1:size(E_total(k,:),2),
    if It(1,i) > a0,
        a0=abs(E_total(k,i)).^2;
    end
end
a1=0.5*a0;
p=1;
for i=1:size(E_total(k,:),2)-1,
    if (It(1,i)-a1)*(It(1,i+1)-a1)<0,
        a3(1,p)=t(1,i);
        a5(1,p)=i;
        p=p+1;
    end
end
a4(k)=a3(1,3)-a3(1,2);
clear a3;
samp=size(t,2);
Spectrum=fft(E_total(k,:),samp);
f=1/samp_time*(1:samp)/samp;
I_f=Spectrum.*conj(Spectrum)/samp;
I_fphase=atan(imag(Spectrum)./real(Spectrum));
\end{verbatim}

A little modification is needed when you want to generate data like FIG. 5. We only change the parameter {\verb amp } and neglect plotting procedure.
  
\end{document}
