# latex_exampleold

preamble: \usepackage{minted}
text:
  To include a python code example, use the following:
  ```
  \begin{minted}{python}
  import numpy as np
  import matplotlib.pyplot as plt

  fig, ax = plt.subplots()
  x = np.linspace(0,10,11)
  ax.plot(x, np.sqrt(x))
  plt.show()
  \end{minted}
  ```


text:
  ```
  \begin{align*}
    E &= mc^2 &\qquad G &= oe^4 \\
    F &= nd^3 &\qquad H &= pf^5
  \end{align*}
  ```


text:
  ```
  \begin{tabular}{lcr}
    \toprule
    left & center & right \\
    \midrule
    1 & 2 & 3 \\
    4 & 5 & 6 \\
    \bottomrule
  \end{tabular}
  ```


text:
  ```
  \captionsetup{
    justification=raggedright,
    singlelinecheck=false,
    margin=10pt,
    font=small,
    labelfont=bf
  }
  ```


text:
  Suppose that you have a figure labeled with the command `\label{fig:gamma}`
  In "normal" LaTeX, you create a cross reference to this figure with a command such as
  ```
  Fig.~\ref{fig:gamma}
  ```
  which would output something like **Fig. 3**. When you use the **hyperref** package to turn cross references into active links, only the number 3 would be part of the clickable link.

  The **cleveref** package improves this system by (a) handling the `Fig.~` part for you. If you load **cleveref** with options `nameinlink` and `capitalize`, it will include the `Fig.~` part in the clickable link. If you prefer to have things like equations and figures referenced with lowercase (e.g. `eq. (3)` instead of `Eq. (3)`, don't use the option `capitalize`.

  ## Defining your own reference types

  Suppose you define your own environment with its corresponding counter. For example, a `problem` environment with associated `\newcounter{problem}`. Inside the problem environment, you define the label with a statement like
  ```
  \label[problem]{prob:<thelabel>}
  ```
  Then you can get **cleveref** to work with problems by putting
  ```
  \crefname{problem}{Problem}{Problems}
  ```
  in your document preamble after loading the **cleveref** package.


text:
  ```
  \renewcommand{\footrulewidth}{0.4pt}
  \newcommand*{\ffcolor}{DarkOliveGreen}
  \newcommand*{\ffstyle}[1]{{\scriptsize\textcolor{\ffcolor}{\nouppercase#1}}}
  \renewcommand{\headrule}{{\color{\ffcolor}%
    \hrule width\headwidth height\headrulewidth \vskip-\headrulewidth}}
  \renewcommand{\footrule}{{\color{\ffcolor}%
    \hrule width\headwidth height\footrulewidth \vskip-\footrulewidth}}

  % Set the general page style
  \fancyhf{} % clear all fields
  \fancyfoot[LO,RE]{\ffstyle{Document name}}
  \fancyfoot[LE,RO]{\ffstyle{\theDate}}
  \fancyfoot[CO,CE]{\ffstyle{\thepage}}

  % Set the style for the first page
  \fancypagestyle{first}{\fancyhf{}%
    \fancyfoot[LO,RE]{\ffstyle{Something different}}%
    \fancyfoot[RO,LE]{\ffstyle{\theDate}}%
    \renewcommand\headrulewidth{0pt}

  % To make the headers and footers show up, we need to set the pagestyles
  \thispagestyle{first}
  \pagestyle{fancy}
  }
  ```


text:
  To link to an external source, you can use the command
  ```
  \hyperref{url}{x}{y}{tag text}
  ```
  which will construct a hyperlink to `url/x.y`. If you don't need a separate tag text from the url itself, you can just use
  ```
  \url{url}
  ```


error: []
text:
  ~~~
  \newverbcommand{〈\macro〉}[〈\verbmacro〉]{〈code before〉}{〈code after〉}
  \renewverbcommand{〈\macro〉}[〈\verbmacro〉]{〈code before〉}{〈code after〉} 
  \provideverbcommand{〈\macro〉}[〈\verbmacro〉]{〈code before〉}{〈code after〉}
  ~~~


text:
  Near the top of the document, include

  ```
  \pagenumbering{arabic} % or other style
  ```

  In a footer, use something like

  ```
  \fancyfoot[C]{\thepage of \pageref{LastPage}}
  ```


text:
  The following surrounds the fraction a/b with parentheses that are sized appropriately. You can also use brackets and braces.
  ```
  $\qty( \frac{a}{b} )$
  $\qty[ \frac{1-x^2}{2x} ]$
  $\qty[ \frac{1-x^2}{2x} ]$
  ```

  For the measure in an integral, use `\dd{x}`, which sets the **d** in roman font and adjust spacing before it according to context. Use `\dv{x}` for d/dx, `\dv{f}{x}` for $\df/dx$, and `\dv*{f}{x}` for an inline version. Similar expressions for partial derivatives use `\pdv`.

  ```
  \begin{equation}
    \Gamma[n+1) = \int_{0}^{\infty} x^n e^{-x} \dd{x}
    \label{eq:gamma}
  \end{equation}
  ```
  For vectors, you can use `\vb{F}` or `\vb*{\sigma)` for vectors set in boldface, `\va{F}` or `\va*{\theta}` for an arrow, and `\vu{a}` or `\vu*{\theta}` for unit vectors. There is also full support for Dirac bra-ket notation.


text:
  ```
  \begin{SCfigure}[5][tb]
    \centering
    \includegraphics[width=3in]{2D-0015-stop}
    \caption{The dependence of stopping condition on steps in a two-dimensional random
    walk as a function of the total number of steps taken
    before encountering an overlap of atoms or bonds.}
    \label{fig:stopping}
  \end{SCfigure}
  ```
  The first optional parameter ([5] in the example) describes the space to reserve for the caption. A large value (e.g., 50) reserves the maximum width that is possible.


text:
  ```
  The counterweight's mass was $m = \SI{2.34(2)e3}{kg}$
  ```
  Note: recent versions of **siunitx** have changed the standard macro to `\qty`, which interferes with the macro with the same name in the **physics** package. If you load **physics** before **siunitx**, the **physics** version is retained and you can use `\SI{val}{units}` as before.

  If you need to put something in the numeric part that collides with **siunitx**'s formatter, you can still use `\unit{units}` to set the units part properly. 


  + For angles, use `\ang{1;2;3}` to get 1°2'3''
  + `\num[options]{number}` for numerical values: `\num[color=blue]{1.234e7}`
  + `\unit[<options>]{<unit>}` for just the unit portion. Separate units with periods: `\unit{kg.m/s^2}`
  + \complexnum[<options>]{<number>], as in `\complexnum{1+-2i}`


text:
  The example below shows how to create a figure with two subfigures. Each subfigure is drawn by **tikz**, but you can also use `\includegraphics[options]{filename}`
  ``` latex
  \begin{figure}[tb]
  	\centering
  		\subcaptionbox{a microstate\label{fig:amicrostate}}{
  			\begin{tikzpicture}[scale=1]
  				\def\x{3}
  				\def\y{3}
  				\draw[dashed] (0,0) -- (0,\y);
  				\draw[line width=4pt,color=SlateGray] (-\x,0) rectangle (\x,\y);
  				\foreach \x/\y/\n in {-2/2.5/1, 1/2.3/2, -0.7/1.3/3, -2.5/0.4/4, 1/1.4/5, 2.6/0.8/6, -2.2/2/7} {
  				  \node[draw,circle,inner sep=2pt,fill=yellow!20] at (\x,\y) {\n};
  				}
  			  \end{tikzpicture}%\includegraphics[scale=.75]{microstate}}
  		}
  		\qquad
  		\subcaptionbox{a macrostate\label{fig:amacrostate}}{
  			\begin{tikzpicture}[scale=1]
  				\def\x{3}
  				\def\y{3}
  				\draw[dashed] (0,0) -- (0,\y);
  				\draw[line width=4pt,color=SlateGray] (-\x,0) rectangle (\x,\y);
  				\foreach \x/\y/\n in {-2/2.5/1, 1/2.3/2, -0.7/1.3/3, -2.5/0.4/4, 1/1.4/5, 2.6/0.8/6, -2.2/2/7} {
  				  \node[draw,circle,inner sep=2pt,fill=yellow!20] at (\x,\y) {\phantom{\n}};
  				}
  				\node[anchor=south] at (-\x/3,0.1) {$N-r$};
  				\node[anchor=south] at (\x/3,0.1) {$\vphantom{N}r$};
  			  \end{tikzpicture}%\includegraphics[scale=.75]{macrostate}}
  		}
  	\caption[Microstates and macrostates of a gas]{
  	(a) A microstate of the gas in the box having $r$ molecules on 
  	the right half of the box and $N-r$ molecules on the left. Note that 
  	we obtain a different microstate on permuting any pair of molecules 
  	of which one is on the left and the other is on the right. (b) The 
  	corresponding macrostate. It is characterized solely by having $r = 
  	3$ molecules on the right, and $N-r=4$ molecules on the left. 
  	}
  	\label{fig:boltzbeetleA}
  \end{figure}

  ```
  By default, the letter of the subfigure is tacked onto the figure number (**Fig. 3a**). If you want the subfigure designation to be set off in parentheses (**Fig. 3(a)**), load the package with the following code:
  ```
  \usepackage[labelformat=simple]{sub caption}
  \renewcommand\thesubfigure{(\alph{subfigure})}
  ```


text:
  ```
  % \begin{wrapfigure}[nlines]{<rl>}[overhang into margin]{width or 0pt}
  \begin{wrapfigure}{r}{0pt}
    \includegraphics[width=0.75\textwidth]{file}
    \caption[short version]{Long version of the caption}
    \label{fig:wrapper}
  \end{wrapfigure}
  ```


error: []
text:
  ``` latex
  \ifthenelse{\isempty{#1}}
      {Do something when the first argument is empty}
      {Do something else when it isn't empty}
  ```


# latex_example

text:
  ```
  \begin{align*}
    E &= mc^2 &\qquad G &= oe^4 \\
    F &= nd^3 &\qquad H &= pf^5
  \end{align*}
  ```


  amssymb
  array
  booktabs
  ```
  \begin{tabular}{lcr}
    \toprule
    left & center & right \\
    \midrule
    1 & 2 & 3 \\
    4 & 5 & 6 \\
    \bottomrule
  \end{tabular}
  ```


preamble:
  \usepackage[nameinlink,capitalize]{cleveref}
  \crefname{problem}{Problem}{Problems}

text:
  Suppose that you have a figure labeled with the command `\label{fig:sample}`, such as by using

  ```
  \begin{figure}[t!]
    \centering
    \includegraphics[width=2in]{sample}
    \caption{This is just an example figure.}
    \label{fig:sample}
  \end{figure}
  ```

  In "normal" LaTeX, you create a cross reference to this figure with a command such as
  ``` x
  Fig.~\ref{fig:gamma}
  ```
  which would output something like **Fig. 3**. When you use the **hyperref** package to turn cross references into active links, only the number 3 would be part of the clickable link.

  The **cleveref** package improves this system by (a) handling the `Fig.~` part for you. If you load **cleveref** with options `nameinlink` and `capitalize`, it will include the `Fig.~` part in the clickable link. If you prefer to have things like equations and figures referenced with lowercase (e.g. `eq. (3)` instead of `Eq. (3)`, don't use the option `capitalize`.

  ```
  So, to recap, if you write \cref{fig:sample}, you will get a
  very clever reference to that figure.

  \lipsum[2]
  ```

  ## Defining your own reference types

  Suppose you define your own environment with its corresponding counter. For example, a `problem` environment with associated `\newcounter{problem}`. Inside the problem environment, you define the label with a statement like

  ```
  \label[problem]{prob:<thelabel>}
  ```
  Then you can get **cleveref** to work with problems by putting
  ```
  \crefname{problem}{Problem}{Problems}
  ```
  in your document preamble after loading the **cleveref** package.


text:
  ```
  \renewcommand{\footrulewidth}{0.4pt}
  \renewcommand*{\ffcolor}{DarkOliveGreen}
  \renewcommand*{\ffstyle}[1]{{\scriptsize\textcolor{\ffcolor}{\nouppercase#1}}}
  \renewcommand{\headrule}{{\color{\ffcolor}%
    \hrule width\headwidth height\headrulewidth \vskip-\headrulewidth}}
  \renewcommand{\footrule}{{\color{\ffcolor}%
    \hrule width\headwidth height\footrulewidth \vskip-\footrulewidth}}

  % Set the general page style
  \fancyhf{} % clear all fields
  \fancyfoot[LO,RE]{\ffstyle{Document name}}
  \fancyfoot[LE,RO]{\ffstyle{the Date}}
  \fancyfoot[CO,CE]{\ffstyle{\thepage}}

  % Set the style for the first page
  \fancypagestyle{first}{\fancyhf{}%
    \fancyfoot[LO,RE]{\ffstyle{Something different}}%
    \fancyfoot[RO,LE]{\ffstyle{theDate}}%
    \renewcommand\headrulewidth{0pt}
  }

  % To make the headers and footers show up, we need to set the pagestyles
  \thispagestyle{first}
  \pagestyle{fancy}

  \lipsum[1]
  ```

preamble:
  \usepackage{fancyhdr}
  \usepackage[width=4in,height=4in]{geometry}
  \usepackage{lipsum}

  \renewcommand{\footrulewidth}{0.4pt}
  \newcommand*{\ffcolor}{DarkOliveGreen}
  \newcommand*{\ffstyle}[1]{{\scriptsize\textcolor{\ffcolor}{\nouppercase#1}}}
  \renewcommand{\headrule}{{\color{\ffcolor}%
    \hrule width\headwidth height\headrulewidth \vskip-\headrulewidth}}
  \renewcommand{\footrule}{{\color{\ffcolor}%
    \hrule width\headwidth height\footrulewidth \vskip-\footrulewidth}}

  % Set the general page style
  \fancyhf{} % clear all fields
  \fancyfoot[LO,RE]{\ffstyle{Document name}}
  \fancyfoot[LE,RO]{\ffstyle{the Date}}
  \fancyfoot[CO,CE]{\ffstyle{\thepage}}

  % Set the style for the first page
  \fancypagestyle{first}{\fancyhf{}%
    \fancyfoot[LO,RE]{\ffstyle{Something different}}%
    \fancyfoot[RO,LE]{\ffstyle{theDate}}%
    \renewcommand\headrulewidth{0pt}
  }
  % To make the headers and footers show up, we need to set the pagestyles
  \pagestyle{fancy}

error:
  ! File ended while scanning use of \f@nch@pagestyle.
  ! Emergency stop.
  !  ==> Fatal error occurred, no output PDF file produced!


preamble: \usepackage{minted}
text:
  To include a python code example, use the following:

  ```
  \begin{minted}{python}
  import numpy as np
  import matplotlib.pyplot as plt

  fig, ax = plt.subplots()
  x = np.linspace(0,10,11)
  ax.plot(x, np.sqrt(x))
  plt.show()
  \end{minted}
  ```


text:
  Near the top of the document, include

  ```
  \pagenumbering{arabic} % or other style 
  ```

  In the preamble to adjust a footer, use something like
  `\fancyfoot[C]{\thepage\ of \pageref{LastPage}}`

  ```
  \thispagestyle{fancy}
  \lipsum[1]
  ```

preamble:
  \usepackage{fancyhdr}
  \usepackage{pageslts}
  \pagenumbering{arabic}
  \usepackage[width=4in, height=3.5in]{geometry}
  \usepackage{lipsum}

  \fancyfoot[C]{\thepage\ of \pageref{LastPage}}
  \pagestyle{fancy}


preamble: \usepackage{sidecap}
text:
  ```
  \begin{SCfigure}[5][tb]
    \centering
    \includegraphics[width=2in]{sample}
    \caption{The dependence of stopping condition on steps in a two-dimensional random
    walk as a function of the total number of steps taken
    before encountering an overlap of atoms or bonds.}
    \label{fig:stopping}
  \end{SCfigure}

  \lipsum[1-2]
  ```
  The first optional parameter ([5] in the example) describes the space to reserve for the caption. A large value (e.g., 50) reserves the maximum width that is possible.


preamble: \usepackage{subcaption}
text:
  The example below shows how to create a figure with two subfigures. Each subfigure is drawn by **tikz**, but you can also use `\includegraphics[options]{filename}`

  ```
  \begin{figure}[tb]
  	\centering
  		\subcaptionbox{a microstate\label{fig:amicrostate}}{
  			\begin{tikzpicture}[scale=1]
  				\def\x{3}
  				\def\y{3}
  				\draw[dashed] (0,0) -- (0,\y);
  				\draw[line width=4pt,color=SlateGray] (-\x,0) rectangle (\x,\y);
  				\foreach \x/\y/\n in {-2/2.5/1, 1/2.3/2, -0.7/1.3/3,
                                        -2.5/0.4/4, 1/1.4/5, 2.6/0.8/6, -2.2/2/7} {
  				  \node[draw,circle,inner sep=2pt,fill=yellow!20] at (\x,\y) {\n};
  				}
  			  \end{tikzpicture}
  		}
  		\qquad
  		\subcaptionbox{a macrostate\label{fig:amacrostate}}{
  			\begin{tikzpicture}[scale=1]
  				\def\x{3}
  				\def\y{3}
  				\draw[dashed] (0,0) -- (0,\y);
  				\draw[line width=4pt,color=SlateGray] (-\x,0) rectangle (\x,\y);
  				\foreach \x/\y/\n in {-2/2.5/1, 1/2.3/2, -0.7/1.3/3,
                                       -2.5/0.4/4, 1/1.4/5, 2.6/0.8/6, -2.2/2/7} {
  				  \node[draw,circle,inner sep=2pt,fill=yellow!20] at (\x,\y) {\phantom{\n}};
  				}
  				\node[anchor=south] at (-\x/3,0.1) {$N-r$};
  				\node[anchor=south] at (\x/3,0.1) {$\vphantom{N}r$};
  			  \end{tikzpicture}
  		}
  	\caption[Microstates and macrostates of a gas]{
  	(a) A microstate of the gas in the box having $r$ molecules on 
  	the right half of the box and $N-r$ molecules on the left. Note that 
  	we obtain a different microstate on permuting any pair of molecules 
  	of which one is on the left and the other is on the right. (b) The 
  	corresponding macrostate. It is characterized solely by having $r = 
  	3$ molecules on the right, and $N-r=4$ molecules on the left. 
  	}
  	\label{fig:boltzbeetleA}
  \end{figure}

  ```

  By default, the letter of the subfigure is tacked onto the figure number (**Fig. 3a**). If you want the subfigure designation to be set off in parentheses (**Fig. 3(a)**), load the package with the following code:

  ``` x
  \usepackage[labelformat=simple]{subcaption}
  \renewcommand\thesubfigure{(\alph{subfigure})}
  ```


preamble:
  \usepackage{wrapfig}
  \usepackage{lipsum}

text:
  ```
  % \begin{wrapfigure}[nlines]{<rl>}[overhang into margin]{width or 0pt}
  \begin{wrapfigure}{r}{0pt}
    \includegraphics[width=2in]{sample}
    \caption[short version]{Long version of the caption}
    \label{fig:wrapper}
  \end{wrapfigure}
  \lipsum[1-2]
  ```


preamble: \usepackage[svgnames]{xcolor}
text:
  ```
  There are two ways to change the color of text. To alter the color 
  in a small block, you can use \textcolor{DarkRed}{the textcolor 
  command}. To change the color more generally, use \color{DarkGreen}
  a color command, which stays in effect until you exit the current scope. 
  ```


preamble: NOCOMPILE
text:
  ```
  \definecolor{mygray}{gray}{0.4}
  \definecolor{darkgray}{gray}{0.25}
  \definecolor{gold}{RGB}{178, 139, 51} %hmc-official gold

  \usepackage{titletoc}
  \titlecontents{chapter}[3em]{\addvspace{1.5pc}\bfseries}%
  {\contentslabel{3.em}}{\sf}%
  {\hspace{0.5em}{\color{gold}\footnotesize$\bullet$}\hspace{0.5em}{\thecontentspage}}
  \titlecontents{section}[3em]{}%
  {\contentslabel{3.em}}{}%
  {\hspace{0.5em}{\color{mygray}\footnotesize$\bullet$}\hspace{0.5em}{\thecontentspage}}
  \titlecontents{subsection}[3em + 3.5em]{\small}%
  {\contentslabel{3em}}{ }%
  {\hspace{0.5em}{\color{darkgray}\scriptsize$\bullet$}\hspace{0.5em}{\thecontentspage}}
  \titlecontents{figure}[3em]{}%
  {\contentslabel{3.em}}{\sf}%
  {\hspace{0.5em}{\color{mygray}\footnotesize$\bullet$}\hspace{0.5em}{\thecontentspage}}
  \titlecontents{table}[3em]{}%
  {\contentslabel{3.em}}{\sf}%
  {\hspace{0.5em}{\color{mygray}\footnotesize$\bullet$}\hspace{0.5em}{\thecontentspage}}
  ```


preamble: NOCOMPILE
text:
  ```
  \titleformat{\chapter}[display]
    {\normalfont\bfseries\color{gold}}
    {\filleft\hspace*{-60pt}%
      \rotatebox[origin=c]{90}{%
        \color{black}\Large%
          \normalfont\textsc\chaptertitlename%
      }\hspace{10pt}%
      {\setlength\fboxsep{0pt}%
      \colorbox{gold}{\parbox[c][3cm][c]{2.5cm}{%
        \centering\color{white}\fontsize{80}{90}\sansserif\thechapter}%
      }}%
    }
    {10pt}
    {\color{gold}\titlerule[2pt]\vskip3pt
    \color{black}\titlerule[1pt]\vskip4pt
    \color{gold}\huge\sansserif}


  \titleformat{\section}
    {\sansserif\Large\bfseries\color{mygray}}
    {\thesection}{1em}{}

  \titleformat{\subsection}
    {\color{darkgray}\sansserif\large\bfseries}
    {\thesubsection}{1em}{}

    \titleformat{\subsubsection}
    {\sansserif\bfseries}
    {\thesubsection}{1em}{}{}
  ```


preamble:
  \usepackage{lipsum}
  \usepackage{multicol}
  \usepackage[width=4in,height=7in]{geometry}

text:
  ```
  This part is set before we start using multicols. The text spans the width of the page.
  \small
  \begin{multicols}{2}
  \lipsum[1]
  \end{multicols}

  And now with the star:

  \begin{multicols*}{2}
  \lipsum[2]
  \end{multicols*}
  ```


text:
  In the preamble:
  ``` x
  \usepackage{caption}
  \DeclareCaptionFont{quack}{\sffamily}
  \captionsetup{
    font={quack,small},
    labelfont=bf,
    justification=raggedright,
    singlelinecheck=false,
    margin=10pt
  }
  ```

  In the document:

  ```
  \begin{figure}[b!]
    \centering
    \includegraphics[width=0.5\textwidth]{sample}
    \caption{This is my lovely caption.}
    \label{fig:caption}
  \end{figure}
  \lipsum[1]
  ```

preamble:
  \usepackage{caption}
  \DeclareCaptionFont{quack}{\sansserif}
  \captionsetup{font={quack,small}, labelfont=bf}
  \usepackage[width=4in, height=6in]{geometry}


# mpl_graph

created: 2022-06-10 21:49:32.537448
modified: 2022-06-10 22:23:12.672898
description: Shaded region included in a plot of binomial and normal distribution
mime: webp
code:
  import matplotlib
  plt.rc('text.latex', preamble=r'\usepackage{amsmath}')

  def binomgraph(N:int, p:float, ymax=None):
      """Generate a figure showing a comparison of the
      binomial distribution and the gaussian approximation
      with shading to indicate the region within one standard
      deviation.
      """
      from scipy.stats import binom, norm
      q = 1 - p
      sigma = np.sqrt(N * p * q)
      mean = N * p
      x = np.linspace(0, N, N+1)
      y = binom.pmf(x, N, p)
      
      xnorm = np.linspace(mean-sigma, mean+sigma, 60)
      ynorm = norm.pdf(xnorm, mean, sigma)
      fig, ax = plt.subplots(figsize=(4,3))
      ax.step(x, y, where='mid', color=traces[0])
      ax.fill_between(xnorm, ynorm, color=traces[0], alpha=0.25)
      ax.set_ylim(0)
      ax.set_xlim(0,N)
      ax.set_xlabel("$n$")
      ax.set_ylabel(r"$\mathcal{P}$")
      xlabel, ylabel = 0.1*N, 0.8*np.max(y)
      ax.text(xlabel, ylabel, r"\begin{align*}N &= %d \\ p &= %g \\ \sigma &= %.2f\end{align*}" % (N, p, sigma))
      ax.axvline(N*p, ls="--", c='k')
      fig.subplots_adjust(bottom=0.15, left=0.15)

  binomgraph(25, 0.6)


created: 2022-06-10 23:37:06.360345
modified: 2022-06-10 23:37:06.360456
description: Simple sine
mime: webp
code:
  t = np.arange(0, 10, 0.1)
  fig, ax = plt.subplots()
  omega = 2 * np.pi / 5
  ax.plot(t, np.sin(omega * t))
  ax.plot(t, np.cos(omega * t), 'ro')
  ax.set_xlabel("$t$")
  ax.set_ylabel(r"$\sin(\omega t)$ and $\cos(\omega t)$");


created: 2022-06-10 23:46:52.273293
modified: 2022-06-10 23:46:52.273427
description: Error bars
mime: webp
code:
  import pandas as pd

  the_data = pd.read_table('https://www.physics.hmc.edu/courses/p134/CircularMoore2004.txt')
  fig, ax = plt.subplots()
  ax.errorbar(the_data.x, the_data.I, the_data.I_err*50, fmt='rs', capsize=3, ms=6)
  ax.set_xlabel("$x$ (mm)")
  ax.set_ylabel("$I$ (V)")


created: 2022-07-21 14:25:57.402603
modified: 2022-07-21 14:25:57.402637
description: Shows how to put the legend outside the frame
mime: webp
code:
  def PDE0(x, t, D=1, L=1, n_max=20):
      u = np.zeros(len(x))
      for n in range(n_max):
          m = 2 * n + 1
          u += 4/(np.pi * m) * np.sin(m*np.pi*x/L) * np.exp(-(m*np.pi/L)**2 * D * t)
      return u

  fig, ax = plt.subplots()
  x = np.linspace(0, 1, 201)
  for t in (0, 0.01, 0.02, 0.05, 0.1, 0.2):
      ax.plot(x, PDE0(x, t, n_max=500), label=f"$t = {t}$")
  plt.subplots_adjust(right=0.8)
  ax.legend(loc='center left', bbox_to_anchor=(1,0.5))
  ax.set_xlabel("$x/L$")
  ax.set_ylabel("$u$")


created: 2022-07-23 13:12:39.290156
modified: 2022-07-23 13:12:39.290189
description: Solution to Laplace's equation in a circle
keywords: polar contour
mime: webp
code:
  azimuths = np.radians(np.arange(0, 360.1, 2))
  zeniths = np.linspace(0, 1, 101)

  r, theta = np.meshgrid(zeniths, azimuths)
  V = np.zeros((len(azimuths), len(zeniths)))

  for n in range(3, 99, 6):
      V += 3/n * np.tensordot(
          np.sin(n*azimuths), np.power(zeniths, n), axes=0)

  fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))
  im = ax.contourf(theta, r, V, cmap='gist_ncar', levels=32)
  plt.colorbar(im);


created: 2022-07-23 17:56:19.638442
modified: 2022-07-23 17:56:19.638475
description: Heat map solution to Laplace's equation
mime: webp
code:
  def Lapl2d(nmax):
      x = np.linspace(0, 1, 101)
      y = np.linspace(0, 1, 101)
      V = np.zeros((len(y), len(x)))
      for n in range(1, nmax+1, 2):
          npi = n * np.pi
          c = 4 / (npi * np.sinh(npi))
          a, b = np.sin(npi * x), np.sinh(npi*y)
          V += c*np.tensordot(b, a, axes=0)

      fig = plt.figure()
      im = plt.imshow(V, cmap='gist_ncar', origin='lower', extent=(0,1,0,1),vmin=0, vmax=1.2)
      fig.colorbar(im, label="$V$")
      fig.axes[0].set_title(r"$n_{\rm max} = %d$" % nmax)

  Lapl2d(19)


created: 2022-08-22 00:59:00.900043
modified: 2022-08-22 00:59:43.649282
description: A power-law plot including transparency and overlaid fits.
mime: webp
code:
  from fit import Power
  from math import factorial
  n = np.linspace(1, 50, 50, dtype=np.int32)
  factn = [factorial(x) for x in n]
  est1 = np.power(n / np.math.exp(1), n)
  est2 = est1 * np.sqrt(2 * np.pi * n)
  err1 = 1. - est1/factn
  err2 = 1. - est2/factn
  err3 = 1. - est2 * (1. + 1/(12*n)) / factn
  fig, ax = plt.subplots(figsize=(6,3.5))
  ax.loglog(n, err1, '.', label=r"$n^n e^{-n}$")
  ax.loglog(n, err2, 'o', label=r"$\times\sqrt{2\pi n}$", alpha=0.5)
  ax.loglog(n, err3, 'o', label=r"$\times\left(1 + \frac{1}{12 n}\right)$", alpha=0.5)
  n0 = 20
  pf2 = Power(n[n0:], err2[n0:], p0=(-1, 0.1, 0), hold="001")
  pf3 = Power(n[n0:], err3[n0:], p0=(-1, 1, 0), hold="001")
  ax.loglog(n[n0:], pf2(n[n0:]), color='#A00000', lw=2.5, label=r"$\nu=-0.999$")
  ax.loglog(n[n0:], pf3(n[n0:]), color='#006000', lw=2.5, label=r"$\nu=-1.97$")
  ax.set_xlabel("$n$")
  ax.tick_params(axis='x', pad=8)
  ax.set_ylabel("Relative Error")
  plt.subplots_adjust(right=0.7, bottom=0.15)
  ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))


created: 2022-09-09 03:45:55.164243
modified: 2022-09-09 03:45:55.164280
description: Run Ehrenfest's model simulation and compare to analytic expression
keywords: Ehrenfest analytic
mime: webp
code:
  from numpy.random import default_rng
  rng = default_rng()

  N = 100
  def run(N:int, rounds:int):
      """
      I'll use an array of 100 integers. If the integer is 1,
      the ball is in urn 1. If it is 0, it is in the other urn.
      At each turn we'll pick a random index and flip its identity.
      """
      balls = np.ones(N)
      balls[rng.permutation(N)[:25]] = 0 # set 25 of them to 0
      state = [np.sum(balls),]
      for n in range(rounds):
          i = rng.integers(0, N)
          balls[i] = 1 - balls[i]
          state.append(np.sum(balls))
      return np.array(state)

  runs = np.array([run(N, 100) for x in range(10)])
  av = np.mean(runs, axis=0)
  # Let's also include a prediction from the continuum model
  j = np.arange(0, 100, 1)
  continuum = 50 + 25 * np.exp(-2*j/100)

  fig, ax = plt.subplots()
  for i in range(10):
      ax.plot(runs[i,:], alpha=0.5)
  ax.plot(av, 'k', label='Mean')
  ax.plot(j, continuum, 'r--', label='Continuum')
  ax.legend(loc='upper right')
  ax.set_xlabel('Round')
  ax.set_ylabel('$N$ in Urn 1');


created: 2022-09-16 01:51:34.866734
modified: 2022-09-16 01:51:34.866769
description: A log-log plot with tick labels adjusted.
mime: webp
code:
  R = 13.605693122994 # eV
  kB = 1.380649e-23 / 1.60217663e-19 # eV/K

  def threetosixovertwo(T):
      kT = kB * T
      nvals = np.arange(2,7,1)
      E = -R / np.power(nvals, 2)
      dE = E - E[0]
      P = nvals**2 * np.exp(-dE/kT)
      return (np.sum(P) - P[0]) / P[0]

  fig, ax = plt.subplots()
  T = np.linspace(6000, 12000, 50)
  ratio = np.array([threetosixovertwo(x) for x in T])
  ax.loglog(T, ratio)
  ax.set_xlabel('$T$ (K)')
  ax.set_ylabel(r'Ratio $\displaystyle\frac{P(3 - 6)}{P(2)}$')
  ax.grid(which='both')
  ax.tick_params(axis='x', which='both', pad=8);


created: 2022-09-29 02:29:34.026528
modified: 2022-09-29 02:35:19.230898
description: Rotational partition function of hydrogen
mime: webp
code:
  import pandas as pd

  mc2 = 938e6                    # eV
  hbarc = 197                    # eV nm
  kB = 1.380649e-23 / 1.602e-19  # eV / K

  r = 0.074                      # nm
  m = mc2
  I = m/2 * r**2                 # eV nm^2
  B = hbarc**2 / (2 * I)         # eV
  epsilon0 = B / kB              # K

  def Ztab(T):
      l = np.arange(0, 100, 1)
      τ = kB * T
      df = pd.DataFrame({'ll1': l*(l+1), 'Ω': 2*l+1}, index=l)
      df['E'] = df.ll1 * B
      df['Boltz'] = df.Ω * np.exp(-df.E / τ)
      df['Elevel'] = df['Boltz'] * df.E
      return df

  def ErotH(T):
      df = Ztab(T)
      Epara = np.sum(df.Elevel[0:100:2]) / np.sum(df.Boltz[0:100:2])
      Eortho = np.sum(df.Elevel[1:100:2]) / np.sum(df.Boltz[1:100:2])
      Ecomb = 0.25 * Epara + 0.75 * Eortho
      return Epara, Eortho, Ecomb

  T = np.linspace(20, 400, 200)
  E = np.array([ErotH(x) for x in T])
  dT = 0.01
  EdE = np.array([ErotH(x+dT) for x in T])
  C = (EdE - E)/(kB*dT)

  fig, ax = plt.subplots(figsize=(5,3.75))
  for n in range(3):
      ax.plot(T, C[:,n], c=traces[n])
  ax.set_xlabel("$T$ (K)")
  ax.set_ylabel(r"$C_{\rm rot}/k_{\rm B}$")
  ax.annotate("parahydrogen", (T[120],C[120,0]), 
              color=traces[0],
              xytext=(15, 15), textcoords='offset points',
              arrowprops=dict(arrowstyle="->",color=traces[0]))
  ax.annotate("orthohydrogen", (T[100], C[100,1]),
              xytext=(20,-40), textcoords='offset points',
              color=traces[1],
              arrowprops=dict(arrowstyle="->", color=traces[1]))
  ax.annotate("combined", xy=(T[130], C[130,2]),
              xytext=(-80,20), textcoords='offset points',
              color=traces[2],
              arrowprops=dict(arrowstyle="->", color=traces[2],
                             ))
  ax.set_xlim(0,400)
  ax.set_ylim(0,1.5)


created: 2022-09-30 16:24:11.794142
modified: 2022-09-30 16:40:22.372723
description: Anharmonic oscillator
keywords: align y axis labels
mime: webp
code:
  def Eavg(x, n_max=15):
      """Compute the average energy in an anharmonic oscillator, in multiples of ε,
      for x = kT / ε."""
      n = np.arange(0, n_max+0.5, 1)
      e = 1.03 * n - 0.03 * n**2
      boltz = np.exp(-e/x)
      Z = np.sum(boltz)
      return Z, np.sum(boltz * e) / Z

  def Cavg(x, n_max=15):
      dx = max(1e-4*x, 1e-6)
      Z, E = Eavg(x, n_max)
      Z1, E1 = Eavg(x+dx, n_max)
      return Z, E, (E1-E)/dx

  xx = np.linspace(0, 3, 50)
  vals5 = np.array([Cavg(y, 5) for y in xx]) # Z, E, C
  vals15 = np.array([Cavg(y, 15) for y in xx]) # Z, E, C
  valsSHO = np.array([[1/(1 - np.exp(-1/x)), 1. / (np.exp(1/x)-1), (0.5/x/np.sinh(0.5/x))**2] for x in xx])
  valsSHO[0,0] = 1.0
  valsSHO[0,1:] = 0.0
  fig, axs = plt.subplots(nrows=3, sharex=True, figsize=(5,5))
  for n in range(3):
      axs[n].plot(xx, valsSHO[:,n], c=traces[0])
      axs[n].plot(xx, vals5[:,n], '-.', c=traces[1])
      axs[n].plot(xx, vals15[:,n], '--', c=traces[2])
  axs[2].set_xlabel(r'$k_{\rm B}T / \epsilon$')
  axs[0].set_ylabel('$Z$')
  axs[1].set_ylabel('$E / \hbar \omega$')
  axs[2].set_ylabel(r'$C/k_{\rm B}$')
  axs[0].set_ylim(1,4)
  axs[1].set_ylim(0,3)
  axs[2].set_ylim(0, 1.5)
  axs[2].set_xlim(0,3)
  axs[2].tick_params(axis='x', which='both', pad=8)
  fig.align_ylabels(axs)


created: 2022-10-01 20:49:29.799804
modified: 2022-10-01 20:49:29.799832
description: Bar chart
mime: webp
code:
  def barplot(ax, x, y, color=traces[0]):
      # modify color to fade interior color
      c = np.array(color) * 0.25 + 0.75
      ax.bar(x, y, width=0.99, align='edge', color=c, edgecolor=color)
      ax.set_xlim(0)
      #ax.fill_between(x, y, color=color, alpha=0.25)

  fig, axs = plt.subplots(ncols=2, figsize=(8.5,3))

  x = 0.25
  l = np.arange(0, 8, 1)
  lfine = np.arange(0, 8, 0.1)
  z = (2 * l + 1) * np.exp(-l*(l+1)*x)
  zfine = (2 * lfine + 1) * np.exp(-lfine*(lfine+1)*x)
  ax = axs[0]
  barplot(ax, l, z)
  ax.plot(lfine, zfine, color=traces[1])
  ax.set_xlabel('$l$')
  ax.set_ylabel(r'$(2l+1) e^{-B l(l+1)}/k_{\rm B}T$')
  ax.set_xlim(0,8)
  ax.set_ylim(0,2)

  x = 0.0015
  l = np.arange(0, 80, 1)
  lfine = np.arange(0, 80, x)
  z = (2 * l + 1) * np.exp(-l*(l+1)*x)
  zfine = (2 * lfine + 1) * np.exp(-lfine*(lfine+1)*x)
  ax = axs[1]
  barplot(ax, l, z)
  ax.plot(lfine, zfine, color=traces[1])
  ax.set_xlabel('$l$')
  #ax.set_ylabel(r'$(2l+1) e^{-\beta B l(l+1)}$')
  ax.set_xlim(0, 80)
  plt.subplots_adjust(bottom=0.15)
  axs[0].text(4, 1.5, r'$B > k_{\rm B}T$')
  axs[0].text(7.8, 1.9, '(a)', horizontalalignment='right',
              verticalalignment='top')
  axs[1].text(45, 17, r'$B \ll k_{\rm B}T$')
  axs[1].text(78, 22, '(b)', horizontalalignment='right',
              verticalalignment='top')


created: 2022-10-02 00:45:04.870620
modified: 2022-10-02 00:46:01.497355
description: HD rotation
mime: webp
code:
  x = 0.074         # nm
  mp = 938e6        # eV/c^2
  mu = 2 * mp / 3
  hbarc = 197       # eV nm
  B = hbarc**2 / (2 * mu * x**2) # eV
  kB = 8.617e-5     # eV / K
  BK = B / kB       # K

  def EHD(T, n_max=8):
      x = BK/T
      l = np.array(range(n_max))
      e = l * (l+1) * x
      z = np.exp(-e) * (2*l+1)
      return T * np.sum(e * z) / np.sum(z)

  def CHD(T):
      dT = 1e-4 if T == 0 else 1e-4 * T
      e1, e2 = EHD(T), EHD(T+dT)
      return (e2 - e1) / dT

  T = np.exp(np.linspace(np.log(20), np.log(400), 81))
  from matplotlib.ticker import ScalarFormatter
  fig, ax = plt.subplots(figsize=(5, 3.75))
  ax.plot(T, [CHD(x) for x in T])
  ax.set_xscale('log')
  ax.set_xlabel('$T$ (K)')
  ax.set_ylabel(r'$C / C_{\rm cl}$')
  ax.set_xticks([20, 50, 100, 200, 400])
  ax.get_xaxis().set_major_formatter(ScalarFormatter())
  ax.set_ylim(0, 1.2)


created: 2022-10-02 16:03:11.418603
modified: 2022-10-02 16:07:27.848943
description: Hydrogen rotational heat capacity
mime: webp
code:
  mc2 = 938e6                    # eV
  hbarc = 197                    # eV nm
  kB = 1.380649e-23 / 1.602e-19  # eV / K

  r = 0.074                      # nm
  m = mc2
  I = m/2 * r**2                 # eV nm^2
  B = hbarc**2 / (2 * I)         # eV
  epsilon0 = B / kB              # K

  def Ztab(T):
      l = np.arange(0, 100, 1)
      τ = kB * T
      df = pd.DataFrame({'ll1': l*(l+1), 'Ω': 2*l+1}, index=l)
      df['E'] = df.ll1 * B
      df['Boltz'] = df.Ω * np.exp(-df.E / τ)
      df['Elevel'] = df['Boltz'] * df.E
      return df

  def ErotH(T):
      df = Ztab(T)
      Epara = np.sum(df.Elevel[0:100:2]) / np.sum(df.Boltz[0:100:2])
      Eortho = np.sum(df.Elevel[1:100:2]) / np.sum(df.Boltz[1:100:2])
      Ecomb = 0.25 * Epara + 0.75 * Eortho
      Efast = np.sum(df.Elevel) / np.sum(df.Boltz)
      return Epara, Eortho, Ecomb, Efast

  T = np.linspace(20, 400, 200)
  E = np.array([ErotH(x) for x in T]) # Epara, Eortho, Ecomb, Efast
  dT = 0.01
  EdE = np.array([ErotH(x+dT) for x in T])
  C = (EdE - E)/(kB*dT)

  fig, ax = plt.subplots(figsize=figsize)
  for n in range(4):
      ax.plot(T, C[:,n], c=traces[n])
  ax.set_xlabel("$T$ (K)")
  ax.set_ylabel(r"$C_{\rm rot}/k_{\rm B}$")
  ax.annotate("parahydrogen", (T[120],C[120,0]), 
              color=traces[0],
              xytext=(15, 15), textcoords='offset points',
              arrowprops=dict(arrowstyle="->",color=traces[0]))
  ax.annotate("orthohydrogen", (T[100], C[100,1]),
              xytext=(20,-40), textcoords='offset points',
              color=traces[1],
              arrowprops=dict(arrowstyle="->", color=traces[1]))
  ax.annotate("combined", xy=(T[90], C[90,2]),
              xytext=(-40,20), textcoords='offset points',
              color=traces[2],
              arrowprops=dict(arrowstyle="->", color=traces[2]))
  ax.annotate("fast\nconversion", xy=(T[22], C[22,3]),
              xytext=(-25,30), textcoords='offset points',
              color=traces[3],
              arrowprops=dict(arrowstyle="->", color=traces[3]))
  ax.set_xlim(0,400)
  ax.set_ylim(0,1.5)


created: 2022-12-13 05:10:57.325918
modified: 2022-12-13 18:13:25.655489
description: Curves of heat capacity as a function of temperature for several simple systems.
mime: webp
code:
  from scipy.integrate import quad, solve_ivp
  def kern(t, y):
      return 0 if t == 0 else t**3 / (np.exp(t)-1)
  sol = solve_ivp(kern, (0,15), [0], dense_output=True)

  mlimit = 6 * np.pi**4/90
  def energy(T):
      x = 1/T
      m = mlimit if x > 15 else sol.sol(x)[0]
      return m * T**4

  def C(T):
      eps = 1e-4
      if isinstance(T, (float,int)):
          return 3*(energy(T+eps)-energy(T))/eps
      else:
          return np.array([3*(energy(x+eps)-energy(x))/eps for x in T])

  temp = np.concatenate((np.arange(0.0,0.1, 0.0025), np.arange(0.1,1.5,0.01)))
  ltemp = np.power(10.0, np.linspace(-1.5,0.3,100))
  curves = [
      lambda x: (2*x*np.sinh(0.5/x))**-2,
      lambda x: (2*x*np.cosh(0.5/x))**-2,
      lambda x: C(x),
      lambda x: x**3
  ]
  fig, ax = plt.subplots(1,2,sharex=False, figsize=(8,4))
  ax[1].set_ylim(1e-3,2)
  ax[0].set_ylim(0,1.2)
  ax[0].set_xlim(0,1.5)
  ax[0].set_xlabel(r'$T/T_0$')
  ax[0].set_ylabel('$C$')
  ax[1].set_xlabel(r'$T/T_0$')
  ls = ['solid', 'dotted','dashed', 'dashdot']
  for n, f in enumerate(curves):
      lab = f"({chr(ord('a')+n)})"
      ax[0].plot(temp, f(temp), ls=ls[n], label=lab)
      ax[1].loglog(ltemp, f(ltemp), ls=ls[n], label=lab)
  ax[1].legend(loc='lower right')
  ax[0].legend(loc='center right');


created: 2022-12-13 05:18:11.271087
modified: 2022-12-13 05:18:11.271116
description: van der Waals curves
mime: webp
code:
  def p(t, v):
      """Compute the van der Waals pressure in reduced variables t and v"""
      return (8*t/3)/(v-1/3) - 3/(v*v)

  plt.close('all')
  fig, ax = plt.subplots()
  ax.set_ylim(-0.5, 3)
  ax.set_xlim(0.4,3.5)
  vols = np.linspace(0.35, 3.5, 121)
  for t in np.arange(1.1, 0.79, -0.05):
      ax.plot(vols, p(t, vols), label=f"{t:.2f}")
  ax.set_xlabel(r"$\tilde{v}$")
  ax.set_ylabel(r"$\tilde{p}$")
  ax.legend();


created: 2022-12-13 05:19:26.659239
modified: 2022-12-13 05:19:26.659272
description: Maxwell construction of coexistence curves
mime: webp
code:
  from scipy.optimize import root_scalar

  class Coexistence:
      """Given a reduced temperature t, find the pressure and the two specific volumes that yield
      the same chemical potential. The method optimize() solves for the coexistence pressure
      by finding the local pressure minimum (p_min) and maximum (p_max) on the isotherm, and
      then using a binary search between these values to find where the chemical potential
      of the liquid is equal to the chemical potential of the vapor. At each guessed coexistence
      pressure, p_target, the volume of the liquid (v_liq) and the volume of the vapor 
      (v_vap) are found by executing a crude search to bracket the root, followed by a 
      call to scipy.optimize.root_scalar, which uses Brent's method. The Gibbs free energy
      per particle is then computed at each of these volumes and the difference in values
      is used to which bound (p_lo or p_hi) to replace with p_target, which is then recomputed
      as the average of the bounds. Iteration stops when the difference in chemical potential
      values has a magnitude smaller than error.
      
      Internal values:
          v_one           volume at which pressure = 1
          v_peak          volume of the local pressure peak
          v_liq           volume (per particle) of the liquid phase
          v_vap           volume (per particle) of the vapor phase
          p_min           pressure minimum (or zero, if it would be negative)
          p_max           pressure at the local pressure peak
          p_lo            initialized to p_min, it becomes the moving lower bound
          p_hi            initialized to p_max, it becomes the moving upper bound
          p_target        always updated to be 0.5 * (p_lo + p_hi)
          
          mu_liquid       Gibbs per particle of the liquid phase
          mu_vapor        Gibbs per particle of the vapor phase
      
      Callable methods:
          __str__         a description of the coexistence values
          curve           arrays of volume and corresponding pressure values along the isotherm
          maxwell         a plot of the Maxwell construction
      """
      
      
      def __init__(self, t, error=1e-6):
          assert t > 0 and t < 1.0
          self.t = t
          self.error = error
          self.v_one = None
          self.p_min = None
          self.v_min = None
          self.p_max = None
          self.v_peak = None
          self.v_liq = None
          self.v_vap = None
          self.set_pressures()
          self.p_lo = self.p_min
          self.p_hi = self.p_max
          self.p_target = 0.5 * (self.p_min + self.p_max)
          self.f = lambda v: self.p(v) - self.p_target
          self.optimize()
      
      def __str__(self):
          return f"For T = {self.t:.3f}, p = {self.p_target:.4f}, v_liq = {self.v_liq:.3f}, v_vap = {self.v_vap:.3f}, mu = {self.mu_vapor:.3f}"
      
      def curve(self):
          """
          Compute the pressure as a function of v
          """
          vvals, pvals = [self.v_liq], [self.p_target]
          # Start at the liquid point and go until we hit p = 2
          dv = -0.01
          while pvals[-1] < 2:
              vvals.append(vvals[-1] + dv)
              pvals.append(self.p(vvals[-1]))
          # now jump across the coexistence region
          vapvals = np.linspace(self.v_vap, max(5, self.v_vap+3), 30)
          vvals.extend(vapvals)
          pvals.extend(self.p(vapvals))
          vvals = np.array(vvals)
          pvals = np.array(pvals)
          n = np.argsort(vvals)
          return vvals[n], pvals[n]

      def maxwell(self):
          """Compute a graph showing the Maxwell construction"""
          fig, ax = plt.subplots()
          ax.set_xlabel(r"$\tilde{v}$")
          ax.set_ylabel(r"$\tilde{p}$")
          vmid = 0.85
          v = np.union1d(np.linspace(self.v_one, vmid, 50),
                         np.linspace(vmid, self.v_vap + 0.5, 50))
          p = self.p(v)
          ax.plot(v, p, 'k') # the isotherm
          # Find the value of v where the pressure crosses
          res = root_scalar(
              self.f,
              bracket=(self.v_min + 0.2*(self.v_peak-self.v_min), self.v_peak))
          v_cross = res.root
          v = np.linspace(self.v_liq, v_cross, 40)
          p = self.p(v)
          p0 = self.p_target * np.ones(len(v))
          ax.fill_between(v, p, p0, color="r", alpha=0.2)
          v = np.linspace(v_cross, self.v_vap, 40)
          p = self.p(v)
          ax.fill_between(v, p0, p, color='b', alpha=0.2)

      
      def optimize(self):
          while True:
              self.v_liquid()
              self.v_vapor()
              dmu = self.mu_liquid - self.mu_vapor
              if abs(dmu) < self.error:
                  return
              if dmu > 0:
                  self.p_lo = self.p_target
              else:
                  self.p_hi = self.p_target
              self.p_target = 0.5 * (self.p_lo + self.p_hi)
          
      def v_liquid(self):
          """Compute the volume of the liquid phase at p_target"""
          vliq = np.linspace(self.v_one, self.v_min, 20)
          pliq = self.f(vliq)
          pindex = np.nonzero(pliq > 0)[0][-1]
          vlo, vhi = vliq[pindex], vliq[pindex+1]
          res = root_scalar(self.f, bracket=(vlo, vhi))
          self.v_liq = res.root
          self.mu_liquid = self.gibbs(self.v_liq)
      
      def v_vapor(self):
          p0, v = self.f(self.v_peak), self.v_peak
          n, dv = 0, 0.25
          while n < 1000:
              v += dv
              n += 1
              p1 = self.f(v)
              if p1 < 0:
                  break
              if n % 10 == 0:
                  dv *= 4            
                  
          res = root_scalar(self.f, bracket=(v-dv, v))
          self.v_vap = res.root
          self.mu_vapor = self.gibbs(self.v_vap)
          
          
      def set_pressures(self):
          # Find a rough maximum pressure for v > 1
          v = 1
          dv = 0.1 if self.t < 0.98 else 0.1 * (1-t)
          p0 = self.p(v)
          while True:
              v += dv
              p1 = self.p(v)
              if p1 < p0:
                  break
              p0 = p1
          self.p_max = p0
          self.v_peak = v - 0.5*dv
          
          # find a rough minimum pressure

          v = 1
          p0 = self.p(v)
          while True:
              v -= 0.01
              p1 = self.p(v)
              if p1 > p0:
                  if self.p_min is None:
                      self.p_min = max(0, p0)
                      self.v_min = v
                  elif p1 >= 1:
                      self.v_one = v # the minimum volume we need consider
                      return
              p0 = p1
      
      def p(self, v):
          """Compute the van der Waals pressure in reduced variables t and v"""
          return (8*self.t/3)/(v-1/3) - 3/(v*v)
      
      def gibbs(self, v):
          """Compute the Gibbs free energy per particle as a function of reduced temperature and volume"""
          return self.t * (v / (v-1/3) - np.log(v-1/3)) - 9 / (4 * v)

  c = Coexistence(0.95)
  c.maxwell()


created: 2022-12-13 05:20:41.407628
modified: 2022-12-13 05:20:41.407660
description: Another approach to coexistence curves
mime: webp
code:
  from scipy.optimize import root_scalar

  class Coexistence:
      """Given a reduced temperature t, find the pressure and the two specific volumes that yield
      the same chemical potential. The method optimize() solves for the coexistence pressure
      by finding the local pressure minimum (p_min) and maximum (p_max) on the isotherm, and
      then using a binary search between these values to find where the chemical potential
      of the liquid is equal to the chemical potential of the vapor. At each guessed coexistence
      pressure, p_target, the volume of the liquid (v_liq) and the volume of the vapor 
      (v_vap) are found by executing a crude search to bracket the root, followed by a 
      call to scipy.optimize.root_scalar, which uses Brent's method. The Gibbs free energy
      per particle is then computed at each of these volumes and the difference in values
      is used to which bound (p_lo or p_hi) to replace with p_target, which is then recomputed
      as the average of the bounds. Iteration stops when the difference in chemical potential
      values has a magnitude smaller than error.
      
      Internal values:
          v_one           volume at which pressure = 1
          v_peak          volume of the local pressure peak
          v_liq           volume (per particle) of the liquid phase
          v_vap           volume (per particle) of the vapor phase
          p_min           pressure minimum (or zero, if it would be negative)
          p_max           pressure at the local pressure peak
          p_lo            initialized to p_min, it becomes the moving lower bound
          p_hi            initialized to p_max, it becomes the moving upper bound
          p_target        always updated to be 0.5 * (p_lo + p_hi)
          
          mu_liquid       Gibbs per particle of the liquid phase
          mu_vapor        Gibbs per particle of the vapor phase
      
      Callable methods:
          __str__         a description of the coexistence values
          curve           arrays of volume and corresponding pressure values along the isotherm
          maxwell         a plot of the Maxwell construction
      """
      
      
      def __init__(self, t, error=1e-6):
          assert t > 0 and t < 1.0
          self.t = t
          self.error = error
          self.v_one = None
          self.p_min = None
          self.v_min = None
          self.p_max = None
          self.v_peak = None
          self.v_liq = None
          self.v_vap = None
          self.set_pressures()
          self.p_lo = self.p_min
          self.p_hi = self.p_max
          self.p_target = 0.5 * (self.p_min + self.p_max)
          self.f = lambda v: self.p(v) - self.p_target
          self.optimize()
      
      def __str__(self):
          return f"For T = {self.t:.3f}, p = {self.p_target:.4f}, v_liq = {self.v_liq:.3f}, v_vap = {self.v_vap:.3f}, mu = {self.mu_vapor:.3f}"
      
      def curve(self):
          """
          Compute the pressure as a function of v
          """
          vvals, pvals = [self.v_liq], [self.p_target]
          # Start at the liquid point and go until we hit p = 2
          dv = -0.01
          while pvals[-1] < 2:
              vvals.append(vvals[-1] + dv)
              pvals.append(self.p(vvals[-1]))
          # now jump across the coexistence region
          vapvals = np.linspace(self.v_vap, max(5, self.v_vap+3), 30)
          vvals.extend(vapvals)
          pvals.extend(self.p(vapvals))
          vvals = np.array(vvals)
          pvals = np.array(pvals)
          n = np.argsort(vvals)
          return vvals[n], pvals[n]

      def maxwell(self):
          """Compute a graph showing the Maxwell construction"""
          fig, ax = plt.subplots()
          ax.set_xlabel(r"$\tilde{v}$")
          ax.set_ylabel(r"$\tilde{p}$")
          vmid = 0.85
          v = np.union1d(np.linspace(self.v_one, vmid, 50),
                         np.linspace(vmid, self.v_vap + 0.5, 50))
          p = self.p(v)
          ax.plot(v, p, 'k') # the isotherm
          # Find the value of v where the pressure crosses
          res = root_scalar(
              self.f,
              bracket=(self.v_min + 0.2*(self.v_peak-self.v_min), self.v_peak))
          v_cross = res.root
          v = np.linspace(self.v_liq, v_cross, 40)
          p = self.p(v)
          p0 = self.p_target * np.ones(len(v))
          ax.fill_between(v, p, p0, color="r", alpha=0.2)
          v = np.linspace(v_cross, self.v_vap, 40)
          p = self.p(v)
          ax.fill_between(v, p0, p, color='b', alpha=0.2)

      
      def optimize(self):
          while True:
              self.v_liquid()
              self.v_vapor()
              dmu = self.mu_liquid - self.mu_vapor
              if abs(dmu) < self.error:
                  return
              if dmu > 0:
                  self.p_lo = self.p_target
              else:
                  self.p_hi = self.p_target
              self.p_target = 0.5 * (self.p_lo + self.p_hi)
          
      def v_liquid(self):
          """Compute the volume of the liquid phase at p_target"""
          vliq = np.linspace(self.v_one, self.v_min, 20)
          pliq = self.f(vliq)
          pindex = np.nonzero(pliq > 0)[0][-1]
          vlo, vhi = vliq[pindex], vliq[pindex+1]
          res = root_scalar(self.f, bracket=(vlo, vhi))
          self.v_liq = res.root
          self.mu_liquid = self.gibbs(self.v_liq)
      
      def v_vapor(self):
          p0, v = self.f(self.v_peak), self.v_peak
          n, dv = 0, 0.25
          while n < 1000:
              v += dv
              n += 1
              p1 = self.f(v)
              if p1 < 0:
                  break
              if n % 10 == 0:
                  dv *= 4            
                  
          res = root_scalar(self.f, bracket=(v-dv, v))
          self.v_vap = res.root
          self.mu_vapor = self.gibbs(self.v_vap)
          
          
      def set_pressures(self):
          # Find a rough maximum pressure for v > 1
          v = 1
          dv = 0.1 if self.t < 0.98 else 0.1 * (1-t)
          p0 = self.p(v)
          while True:
              v += dv
              p1 = self.p(v)
              if p1 < p0:
                  break
              p0 = p1
          self.p_max = p0
          self.v_peak = v - 0.5*dv
          
          # find a rough minimum pressure

          v = 1
          p0 = self.p(v)
          while True:
              v -= 0.01
              p1 = self.p(v)
              if p1 > p0:
                  if self.p_min is None:
                      self.p_min = max(0, p0)
                      self.v_min = v
                  elif p1 >= 1:
                      self.v_one = v # the minimum volume we need consider
                      return
              p0 = p1
      
      def p(self, v):
          """Compute the van der Waals pressure in reduced variables t and v"""
          return (8*self.t/3)/(v-1/3) - 3/(v*v)
      
      def gibbs(self, v):
          """Compute the Gibbs free energy per particle as a function of reduced temperature and volume"""
          return self.t * (v / (v-1/3) - np.log(v-1/3)) - 9 / (4 * v)

  fig, ax = plt.subplots()
  ax.set_ylim(0.15,1.5)
  ax.set_xlabel(r"$\tilde{v}$")
  ax.set_ylabel(r"$\tilde{p}$")
  for t in np.arange(0.799, 0.9991, 0.025):
      c = Coexistence(t)
      vvals, pvals = c.curve()
      ax.plot(vvals, pvals)


created: 2023-02-05 21:32:23.119526
modified: 2023-02-05 21:41:25.763756
description: A sample animation using celluloid.
mime: webp
code:
  from matplotlib.colors import hsv_to_rgb

  from celluloid import Camera

  fig = plt.figure()
  camera = Camera(fig)

  for a in np.linspace(0, 2 * np.pi, 30, endpoint=False):
      x = np.linspace(-3, 3, 800)
      X, Y = np.meshgrid(x, x)
      x = X + 1j * Y
      y = (x ** 2 - 2.5) * (x - 2.5 * 1j) * (x + 2.5 * 1j) \
          * (x - 2 - 1j) ** 2 / ((x - np.exp(1j * a)) ** 2
          * (x - np.exp(1j * 2 * a)) ** 2)

      H = np.angle(y) / (2 * np.pi) + .5
      r = np.log2(1. + np.abs(y))
      S = (1. + np.abs(np.sin(2. * np.pi * r))) / 2.
      V = (1. + np.abs(np.cos(2. * np.pi * r))) / 2.

      rgb = hsv_to_rgb(np.dstack((H, S, V)))
      plt.imshow(rgb)
      camera.snap()
  animation = camera.animate()
  animation


created: 2023-02-05 21:46:59.460765
modified: 2023-03-18 22:49:08.910412
description: Example animation of the convergence of a Fourier series
mime: gif
code:
  from matplotlib.animation import FuncAnimation
  fig, ax = plt.subplots()
  ax.set_xlim(-1, 1)
  ax.set_ylim(-0.2,1.05)
  ax.set_xlabel(r"$t/\tau$")
  ax.set_title("1")
  t = np.linspace(-1,1,201)
  line, = ax.plot(t,np.zeros_like(t))

  def rect(t, n):
      y = np.zeros_like(t)
      y += 1/2
      for i in range(1, n+1):
          bi = (1-(-1)**i) / (np.pi*i)
          y -= bi * np.sin(2 * np.pi * i * t)
      return y

  def frame(n:int):
      y = rect(t, n)
      line.set_ydata(y)
      title = ax.set_title(f"n = {n}")
      return line, title
      
  ani = FuncAnimation(fig, frame, range(50), interval=40, save_count=50, blit=True)
  ani.save("FSconverge.gif", fps=15, dpi=100)


created: 2023-03-06 17:16:27.459619
modified: 2023-03-18 22:02:30.681997
description: Monte Carlo evaluation of pi.
mime: webp
code:
  from numpy.random import default_rng
  rng = default_rng()
  fig, ax = plt.subplots(figsize=(6,6))
  phi = np.linspace(0, np.pi/2, 100)
  ax.plot(np.cos(phi), np.sin(phi), 'k-', lw=2)
  pts = rng.uniform(size=(10000,2))
  r2 = pts[:,0]**2 + pts[:,1]**2
  color = np.where(r2 < 1, 'g', 'r')
  ax.scatter(pts[:,0], pts[:,1], s=3, marker='o', c=color)
  inside = np.count_nonzero(r2<1)
  pi_est = 4 * inside / pts.shape[0]
  ax.set_title(r"$\pi \approx " + "%.4f" % pi_est + "$");


created: 2023-03-18 22:52:58.690841
modified: 2023-03-18 22:56:18.768611
description: Solution to Laplace’s equation in a circular region.
mime: webp
code:
  N = 401
  x = np.linspace(-1.01,1.01,N)
  r, phi = np.zeros((N,N)), np.zeros((N,N))
  for i in range(N):
      x2 = x[i]**2
      for j in range(N):
          r[i,j] = x2 + x[j]**2
          phi[i,j] = np.arctan2(x[i], x[j])
  r = np.sqrt(r)
  phi[r>1] = np.nan
  r[r>1] = np.nan

  V = np.zeros_like(r)
  for n in range(1,49,2):
      V += 4/(n*np.pi) * r**n * np.sin(n*phi)

  plt.close('all')
  im = plt.imshow(V, cmap='gist_ncar', origin='lower', extent=(-1.01,1.01,-1.01,1.01))
  plt.colorbar(im);


created: 2023-03-19 03:12:42.759386
modified: 2023-03-19 03:12:42.759407
description: Electrostatic potential inside and outside a piecewise constant potential along a circle.
mime: webp
code:
  N = 401
  m = 3
  x = np.linspace(-m,m,N)
  r, phi = np.zeros((N,N)), np.zeros((N,N))
  for i in range(N):
      x2 = x[i]**2
      for j in range(N):
          r[i,j] = x2 + x[j]**2
          phi[i,j] = np.arctan2(x[i], x[j])
  r = np.sqrt(r)
  r = np.where(r>1, 1/r, r)

  V = np.zeros_like(r)
  for n in range(1,49,2):
      V += 4/(n*np.pi) * r**n * np.sin(n*phi)

  plt.close('all')
  im = plt.imshow(V, cmap='gist_ncar', origin='lower',
                  extent=(-m,m,-m,m))
  plt.colorbar(im);


created: 2023-03-19 21:34:45.390534
modified: 2023-03-19 21:56:17.844913
description: A random walk
mime: webp
code:
  from numpy.random import default_rng
  rng = default_rng()
  from matplotlib.gridspec import GridSpec

  N = 2000
  loops = 5
  fig = plt.figure(figsize=(6,10))
  gs = GridSpec(loops, 2, width_ratios=[1,2])
  t = np.arange(1, 1+N, 1)

  for loop in range(loops):
      ax = fig.add_subplot(gs[loop,0])
      theta = rng.uniform(-np.pi, np.pi, size=N)
      x, y = np.cumsum(np.cos(theta)), np.cumsum(np.sin(theta))
      r = np.hypot(x,y)
      ax.scatter(x,y, c=t, s=1, cmap='Spectral')
      ax.set_aspect('equal')
      ax = fig.add_subplot(gs[loop,1])
      ax.plot(t, np.sqrt(t), 'k')
      ax.scatter(t,r, c=t, s=1, cmap='Spectral')
      if loop+1 == loops:
          ax.set_xlabel("$n$")


