(* Content-type: application/mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 6.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       145,          7]
NotebookDataLength[     57317,       1564]
NotebookOptionsPosition[     53379,       1441]
NotebookOutlinePosition[     53740,       1457]
CellTagsIndexPosition[     53697,       1454]
WindowFrame->Normal
ContainsDynamic->False*)

(* Beginning of Notebook Content *)
Notebook[{

Cell[CellGroupData[{
Cell["Thermal Equilibrium", "Title"],

Cell["\<\
APAM 1601
Columbia University\
\>", "Subsubtitle"],

Cell[CellGroupData[{

Cell["Introduction", "Section"],

Cell["\<\
In this notebook, we introduce the notion of a \"thermal equilibrium\" and \
investigate simple properties like mean energy and average particle density \
using a biased random walk model.  \
\>", "Text"],

Cell[CellGroupData[{

Cell["A simple atmosphere", "Subsection"],

Cell[TextData[{
 "Our starting point is a model for an ideal atmosphere. The gas within the \
atmosphere is confined, at the top, by a downward pointing gravity and, at \
the bottom, by the surface of the planet. The gas will be characterized by a \
\"temperature\", and the temperature will bias how likely the particle will \
take a random step upward (or downward). If the motion of a gas particle \
causes an increase in the overall energy of the atmosphere, then this motion \
must be restricted by a \"transition probability\" containing the Boltzman \
factor,",
 StyleBox[" Exp[\[Dash]\[CapitalDelta]U/T]",
  FontWeight->"Bold"],
 ", where \[CapitalDelta]",
 StyleBox["U ",
  FontSlant->"Italic"],
 "is the change of energy of the atmosphere."
}], "Text"]
}, Closed]],

Cell[CellGroupData[{

Cell["Monte Carlo Simulation (Metropolis, et al., 1953)", "Subsection"],

Cell[TextData[{
 "We will use the famous \"Metropolis\" algorithm (Metropolis, Rosenbluth, \
Rosenbluth, Teller, and Teller, \"Equation of State Calculations by Fast \
Computing Machines,\" ",
 StyleBox["J. Chem. Phys.",
  FontSlant->"Italic"],
 " ",
 StyleBox["21",
  FontWeight->"Bold"],
 " (1953) pp. 1087-1092) to simulate the thermal fluctuations of our model \
atmosphere. A Monte Carlo simulation uses a random number generator to \
simulate the random fluctuations which occur in a large collection of \
particles."
}], "Text"],

Cell["\<\
The \"Metropolis\" algorithm is one of the is one of the most important and \
widely-used numerical methods in the field of statistical physics.(The story \
goes that the authors searched in vain for another physicist named \
Metropolis,so that the paper could have the symmetrical author list \
Metropolis,Metropolis,Rosenbluth,Rosenbluth,Teller and Teller! Metropolis \
(1915-1999) was an early pioneer in Monte Carlo methods and group leader at \
Los Alamos.His name appears first on the paper (cited above),and as a result \
most people refer to this method as the \"Metropolis\" Monte Carlo \
method.However,Marshall Rosenbluth (1927-2003) actually invented the \
algorithm.Marshall is perhaps one of the greatest pioneers in plasma physics \
and fusion energy science,and most plasma physicists refer to this algorithm \
as the Rosenbluth-Teller-Metropolis method.)\
\>", "Text"],

Cell[GraphicsData["PICT", "\<\
0000000009@0IP0A0_l<0?on0000B00004P000000000U01V0000000N00402P00
0000U01V0380000009@0IP0OP0200800PP0003:X00000@000000000000000000
00400000000000000000040000000000000000000000@0000000U01V00030000
00X02P000000U01V0000OWAYIVH0000000000000001QL71/000000001000IP2D
04P000180000037L004:E4U6AR0XC5YG:@00000000000000000000000000000H
ool0000<LgAbL00006/0000<JFiTN0000040000<M6UVH`0000D00000CDd0:P00
<AB06:1@<Y@F34Z4@V5@PP`dL`n8@f9468A69T4HaTDa]1aeJAmNb6?[AFBE:BM5
bVCYBBZbAb=Rc6@c=Nc5RcB@cJMC5VcfOCJOc5Zd>OD=Z]jTDJUD^SDQ_DfTd>WD
bAb^Fc>XDnVE:WEf_E6Sen`F:`DJOE>UcJcCfVDF]fng`<]G>jGB3`b6`V;a>:gf
<S6nSVoa/Dae1cQNb>Fb^DX]>HnFb>IcK9c:Lc^PCfeDBOIZNIeVg2/eVdF6aEVK
B]BJ_9B7BFCHK7IG3HFjUJ3@icLjJhKKHGBiG>3@J9AJnGjolWThCU3727_XC2Me
FEbOEZC8b34CS<Mg?/g=]GKIWMjGAfnW>KeN_ejG[mNcf3h]kfekjoCf>KdoWkk:
]:fYci_<[SoZnfcA@8PL5@FPb8>JiR<Q12D81S2@@<6`R=^Ja[b9/jQ:<Jai>>^e
ZJ^m4l>Z4gBf[I13GZjnkn:nVK[Y/[jeVJ[knATY4N?`oCk?Docc]50KN@A6l2:=
1Ha;[9bm8J`T;@bSLYb[2d9`dSDZ@i5BU@nZbBYNkBJ<`b[^]jW/Tb?5dfalYbIY
K5k`Z=7LkARo;mb17koCG8do;;=/f:<>M2SE@kPReA:9@/_l/DO:D/bY;CW8fa[a
YS<3[?NDTaUXkK^BlX9RZI?j_:I3dbC?=;`9gCS/9<Tm?96gdNc]?Dkc_?E0@:/M
0ZDd]NE:/53SE9TW2dhJ2``5=VB/5=WdVi[X3g6bT9lVLAA54UEE1<k<=/nMQejK
c:92jlA9GKKEYKJb^Z=E:G^_MC7fX`ejA7MTaC7MeLb5Fdm_kGdnb98dGe<KcP84
^b2XBPbj8C1d8`UAe7fIBL;<<[=AgFmm?Fla4DBmLJ//@a[3G_CKg/Jjj/=dYbK<
FTlAITBU2i^>MZI_JUd<NaRDgkOlnH0f63eoH>Seo75adCQUUHJQ5T>9:<ZH[2l]
`cS:>eZKdJGdUeIdoT4CAATTU:6c33>_M^AA`VKLI3EcEfe^REi@SZ1ifj3ZY?JV
f6;G5Kj9?[ob9`cdEnWfWjV8?6;^9F7^4P^9D[3>;b^45ZHfcaR/EOJG^c/C]g=/
Tj@@l:[/_DVd9fVM?FaeDcfidCYLjSnk8k_Q:FYJR1f;VmlFXTMo?I`SHd5@LT:@
Q>88G:29hLP^WJMbW:J]JT1:AdV`eWeb?i4WCJNdWM<9^TB?YmfO@a5kTBYKS_Ik
XCVK/Ja^LD;9WP9@U=Z9W`;0Dl^4G6dPoj1d45d><A9aQ47V>B2TdiakU29;C>PT
Af1RB?^QA:kEdQi4L=LA:CiCc^11[596ba^ISc6/o4FBeEj[gM;cKhoMW;NbC]iD
:kdPAQXJ8RCZd9GCPcCZ2JhhLlA?BmUoPJ5Yiccg5a?><H4_k6T26UJl[1/0]7I]
aK:hPokgWA4cQhXI@iM6m1k<<]AnY:Ed_k<LciVKM5k4M9JcYo10WOZ5<K7T>I0e
ZAl<J[A=B@H0B5N<dIMlGRS><;g0XPYNgW>>TNAMIQcC3?I6m2:39;BB/]O>Bn;Q
hC<5`C>ChU[=b1T@LN]@^TMXnT2I^R8`dJa5_j4Y6<>Jn6KalSlk`?JaHc@dSe6^
G`NbE`oN:`1PDRCC9[G4fN8hcHXUmHR@n2I6B9A=JS0jJibg<C7@8N/Xd;V__LFl
RN0R1W><a9:T`PLfB44688@F6Q2B1`T<K<DUL>Pa8RUA;1nl_IOZ8Sa@@<DO9oc^
86O=?2NeQAEH6fHJ/SJ:4@<8H4hY3hW>C;jH@357boV=Efn<H[Ji>V/C6HR3Si4F
Y5;0IPULeB:^ARH@PR52PaD986gL@LaHidmUlB]ge1Z4b^I_6Tj5>V5UcTL58a/@
FQ_6L:D^0QAZ:>=TVA<iYOj<PiNM0`_0@E8YHW0=Dnlk9>CVMRjY4l1WTSECRBEV
kSb4c0PL/R<lJWnBnZ@7^HD>CX;5[d8^Wad:n/g;X/EYYLdV4B<8@UX<bDH=2ZTd
XYMGB9Cc2U59[4l25=?8B/`25Xc24[CFCJTdFGB;O>ld]5IKh`8R/eCM1MRE3f53
g;AU;^WMbnQY<F5KOIOFj8kCncA1RoW6^@LYH[odMFDR5IDold`Le=@J@li]bUU5
dZkAl32S`@EUR]9n<5[7c>SYHIf0T7BHT]WNGVKJRR1F3OQ2TCRn;lAa_^gb`L;:
Og4IA<FQ=MP]?FV`e@5942Fd=VG@f@j1:[eO^`LX`UV26^<]2AVkbSB<fV:lZ>e<
H4BgVA?Jhdc/SEfJ^hAMJS=2BdZ8o2nndKA:?b^3;LS^09JF3Wc;jW5=9^g98NH@
PKa4gb4VFV]ikcl379/hiEZeHk_GN]:BM8ULAFFYOJ]f30]5`D]=;9XFQ:jSTMQU
SDak7i?ZNM39c=c?ElRTMe;;79:[QC6KiP2?m1Yi8;bO9DS;oXQH<R1XFQaI;/W8
L[F=2eg[Aj?0QU?BFTe9dRH0_5nB[lf^Oa6JhgDcaZTf]UO<j2:BM9We@CPTKjW?
fYcQ2agD]<KbcSGHM@]VJ2[685X12ZSb4g>VEXJinQiV;ALaY?AnS[Ak8fGAoAnS
J?dGXd@PUIlkD6[A2HoEUJLHdXM5UlgK9B@bdPlD]l^<7ffYcM9anF/8gH0QI[?>
]Q:TZ5^_EhQ5kb<=HDVB^Z6A]Q5TfEY7BNT]Tk>fKU>l2529<:@DHIMa<mg2DeM]
j3BI<Bf_6ZQmUa]jFXhNgU^<4i8/;haXkZ_^]=jIhA1[9n]OERdbW/@:RbE<>0aO
oPkH=TTGl7h;Y3Q?0k_EL8^dm9UPYO>Y5jBfWRn7?KL98io4QgbZBK5I^DdbhJFk
XcIZfnb[giHdagHF_^?=J]gI_QV[NmeY4II__n0:?lSUMfNQ;X71<YlnkZ22S15g
X`?/HP]VlWQNXRakbcCCWmejMK8n9[]K3_bSMH<EDN8^[:Nk1]ZnenNF_e]abVh]
Q_@aXgW:kNkeB7h[23UA:V`dOE@=Q`GPg@N4MhFOT8R]GI]1:l6enOId9JH^RciS
Kdid_6U[K>]eY8CLM?l9lB5o9=fa^YofGdE@7NEmQ9<GNmWRo^<J/GmJW[eoDS;4
QK_G2Nkmgj6L[_Y59[4>`V34PK7O3o1ldai<O5h`^Z?VE6YFRje2</K8nPaTcT_d
aXk4V8]lO^]`ULn/m8n`NJ2D[0LPZb;X;noJLZ9/d8/TbB:lo@j0mTkc16fRDXZf
>DfP^n@T8<G`Yh_j_e0>jNcDlJVPS0O0=`C<eDj/gFlB5Hi>Ki1Pn09BYhc:n`V:
iDV:Z63DPRmB;V;nge0e1>1BFXibo:;8:<`f@U2f10h6oHdU1G352hmZFN8Ta_2B
gVQRDi1lW:h`>n@29jDlA?033XC6hZ58dbDi0F4FW`[f]c0H3fQ`3VI>gVS>UfS:
[^ZH_P2De`34iZDV1B;Y2^?H:bVF:<ildW1Fj4m[4i180a4bXn8P]fl>iDi<oc1V
YBMDjZRd=F_<M2U308eFJni91nONKZ]haa2@n_4:n`m15g44XDj<iXfV2DiR8@FL
@T8T9lkW10VI4lfHoFoBka6W6V@]56MbnY28OZG`cQ5/jTolLhENGR9/QX9G1Zjj
bbDhG`M2eN_diDk70O3Bn`S>`8BH8<F;0/mD[6fT1b9V`J^PVI3R6K18XlXo4fmU
5350mXhE380a3>SHQDk>46a_3e3boc3JDjNn>l<JjDR^<J><jA404FEN96OTohcP
aXYi7S5o5k26[lS8edWLFXcj34/lNO8<0`nn8:d4CjW4@A1G1B1S4dk]1=8FdQ2n
oDhB8TFXeS2;3E5c8]1nlH9/<<A:9F8P_0dX8f]X34<<aXaJoc7ROY42ndSLe/3W
0RULX4F@`8kB[2<2hF8fFXGLJ:VPWE88Y11DDTbU8?8M6S196V<8U/gRiLk2ob=E
3Hdi8f_>A4R^kDiW9/`@i^8PWdc_:JQ^84o08bNLX5:M3j9?2Dn^ULShUHPM46]V
oW4Pi^_4/ZPnYK7jbRf>mU15:A4k4o1HB`oZSPljhVcW<>_ZR`M[5W8d5Y908NW^
_h>Q2U<hgb85<C0@]k96kHi/DT>J8=;8eT`4D;;@V23VJTWY9Rg^>J_0;n>^T@ZZ
ZXY8ma?8_1;o;m2m=`h7=[2k1LcbiM8[5`hXjnDi3V_JcL9V;Y9NP@dD><8BI^[@
[VZbj8228C9bmHFQ6D8JS?9N]`F:/DY`ZF_T34m>oVFN;XXTB<=`GPBm<XoVJ]37
?S=bi^hH8K@I7jX4o_3o1S5c0C8a705k9:8nV:;X[7AA;hDX87@Be^`@8Q>P<JkE
;T1BRF[c0JIa=@ebhLRFmE4P/N8@A4WB;8jdEJlX9j9VZ^LY9c?:8bhLX@S^Yc9R
XC2>YoBJnbgW?a3FRa5^cRjN<<8Q7f8g0fL_2Z>P9VWi:bHY>J2DX59XJP2DehNL
cg9TUa2HF>;W<k<dXhg`2D@/8TI_9jCFa=A9BkBh9]=<2ddI8LkfdN;nU/iE2</0
Mi2?<;>e=6[hkAEBULUX9FM8A4IhDhWh8=4R8f;YB6UO4N8O:28c0/ko@8<5F6oV
hJ;W@aC@IaEb8;FH86RFUDZbb885???DZV:fE6E6<`>_;N8D/LDXZj8PF:Yi5o7W
5niE2?0RQYFIFRl4D:f[3Tfk8^4hIZWhb28b87Dfmo63F01b@YD[ALHV10DVNL;f
<97`Z99S46YX/Jkm>MH08PD^c0//F0F^nD>jhTgOAWBL`6]^gZQ[9QI7;EI5@[7U
4?0MI4Z=IFiH]Fg69BX4/d;fF:KN9jjH8k4@3EHC:hoR[4La=TFKF@^Z8^8PZEFK
Ig9ZmD`]@K6@`>2296R<lLU508<`nZK_Ja8Y:OJfeU=m?_?aK35k5:m70SI=0Tl8
S09GGPX?I23fOH8oDgE//2X=FTY/oVmfZ^@XLY<R]nUm7]DR2TNO2S688ZXFT7Fd
c3L@=aFk09HiAZ]i0E=hnVc@_gLQ<o7S9Lm6YoJe>gI6N4bi>9LUG^9:lT=<@n<K
@c@dZkJFBM2MF><08NYXWh<<^>8:/dYU9f2TZk=KFdEIDg00C8DeLS=i=oL]L_N@
]iK4iF/7G@n_K<m7KPojbkCV_ln^?8Fc8TOZS>j>D>LNJNFA0XYYC;KM;F85<S9_
K^YVNTBLI^Ei87=N=g3/9ON=A[22>[M7NC5cK73naeJi5j^1JhV;1ZM8IZUCMT8:
<`nPO]C?47IK>o4GJL2DJLko9ZObV64FUaIgM`8@`[H8Vm@SH<8fhkLEOT_BA>gA
LY>15?N>_eO[L^a_EW=3COOkJi=50UP24FRi>8<<^Fej10<J:<b`QHUaFHUFQ]5e
Pf@FWh[o;CK;Lg>FebDFTXJaC2BX9KA5LIDj>llIN7QHliO]N;OeQRiLk=Pbm:n^
Q95nmj58l]OX9?EdBX8<Hk>8c9EBX4QXYY5=>E97O<35LmLlQY7_OH2dVjbNFOJT
g<YKLJA@g@c?O_RoA]QKL[AZa_Jo5e3liGP1ELMiGND>b`<@A5AjDY@48jcXIA5k
Tl3D]aF=>b]gO_Km0USGPKRU@dNOG8b6JZ8b?3871fM;LLcFK3M7QGL_TQHjeUE@
V5KHDB]08FbOI]V35J5886>K47LYPOK//oM^cmIo@?GUV:<OIA97K?BPBJ;WAmD4
Vj@/8<AIBhP^eB>fNh>_S=TgUURQI9I1@^aEDNWRZbbG;P;kJJehFIM0jRTiB@8@
@/Yfm7NZ_eG1BZL`;gF88A<aN9TWGCG>47IDIaFI9c@I4Q?<=GQ;DgIPbj_:G>h/
9;VK0_L5WkMf`^FEHAYQYOPUYTZaYYY]U_4H/m0]XToVn4EOP?MPUM:T_/WiF<m`
ZJhL]/3EM7?a7U0URK2HF@o3JQVZ85WKY:gASM7B>_WoG3YWPX_S?0>6NWMMT9O6
PMO1[DFDZJYT8UAdc69?U61B8<[jMVnLaB85PH3WYFFADJ@Fi8iAIeS_[gHmEEEB
XBFA;k8N8^<8M=06MG/Rji5DbjNh9KO2BLDC/g@eDL;[MEW=/o/k]7]3[<FCVg@h
^fBM5Z=FZ4egC9/e>NIE8bF;7S/1KJ3ECDD:;Y>QEIDDULWiPOJAG41C0fLXNaTG
HdFlbfGU3MVU]1^P@GUE]_^S^[]:l3^^;[66356m]J9>F:/mOC4KPoG:D?Oi=0I/
O^Si@_KM46I[Xc61=8D?T5]8g_A=J:gf1R?2O;Y;^HdfRbM3P_EoIg_Y]aRKXgW;
`9C?VGU[@b[`D:QQSi<f1SE_G/FYD1Q4`B;k`5`ISoI487XK>TUaCHX3`K]=4L<3
J>D;Zb>m/XM1<FOJU?IM`?/Ca[aYaY]oRSFKRIa;]a`jUJ>Qa2SfD?^9=gL]>8`8
Vk9c7k[8BIM[F_^m>UWZV;_I9UT6l1]>/g>AU/[2CT=eN`[F9XnPba1Xj]Foa]CK
T1Hn[oPO^1cFV547_G_I`57M>U7cP@ZL9C`R>ZUaDIbi0ZoWPl8Wb_9<GZ/;CG`o
@e/cO2REYMVl3W/QUmaM7BnPjOCmEN`35i3oK>]acCICC?cKdC`mGo`53iFOcc_A
NT>_MC^U_KL3Y/Zkd5>mGaBMLmKIdFJNm`8CF`jcHedUVRnPlFl^j]7WLhk9RO7Y
ESZQIIfII9^4[o2IUE<o^36fQIcgXg_14HoOL7J8;Xn5GmDE`;JA]=YF8CblhmWM
elDcN4K2oi/^Y>59AUS?JkE>m7Le5mW_fGdoa^3W=;e5>Blk_T3Dfb4iS8SOh;8W
SlTST8BIDB>PcQVBX5cU9U[8knZKWKReRf_?dWgJlEgPjOIcgYJcTi_QckgcdaGK
a]T3aY`i_<SGUSLYT]6eI9;AJAIfF;9Z]gPOC=Eo/m[8JM4LNAQ@K7gNn7f7^J?O
NJgW7QQ]g]dgi;_Sf?jS/AFMI5`;ii5k>ba^c?EIUWijX1Eo`j877SKIWcR<RODP
348CR2hdE7Y9V1/]1kaPi9>8fgZO6ggU9IjGSAQaTidkfEfMj]gigk_WPacLiIZM
P4^73oHO^]<T23>mJANk]cHKL3YGFD;[[P=GWJi1RigHU0TogLbfi9g_gUg[JkjE
gai;f=CM61dmWa9Sg3IMa76fnU=hlo/83WlaP/@I]@TT2U/i^]]<ko:jLH96B=JZ
e?1iaSkWaPgDcLgUU@iM9GlEmKK9mEW5mYjYRKigKM;A94l=nVm7XF47n2/G@fJM
Vc0]lK^Ql_YWVID<hja<l[nG^I_mhm`1gQEIm;hB80NcfRh8Ph=0Sg1XDPhA28E3
HQ3hQ4hY43W5c56H_6cW28i5hU5H514G2Xg23E:Hc:c4FYLBYPDYU;YY=IG=IM<R
U<2D>Il@J1@Ie>2e@ic<Yi>Z3=ia>Yl>J2_JTeJXcJ/aJaEZdcJaGF;D[0_EYHk1
Hk=IeIJKEIU9KK?HkKLE9IY9MK]Mk_2994hE98G1[mOl58Y?5XoQhA:hg:CENX9Q
8_4jHF/W:ie?:3=:=A2e@J5<Z3?:O=9iVYeYmAX:1S9I6M<DZ3Ck1E6[F]Y]:m^K
3D[O_EYJmm_[E`5YLDYakkM/5PD7c>M1H?4NS8XCOnY0/?fHoRXo2^>U<70Y3=<]
ZR3?=5?/_<=O[jOXISZIWCMABOTD]I[O96La@:N/cKZZ[J^=b[Zb^2hSQ[B^;R;R
/be@H]:lY6_2j^B^[79:jK`^ZQSa>]3i1XfT;4>`kC2SfkPiY8kk1>nTR=YXf:O/
l8;gaZdZS]PodJ?>V2W_j8;m]JS;lS57;:B:UB</fYbWbN7:^]XfbZBU0AV]f/CP
[H]bhBjgl5[K22i59<IB>Hkl<[g4<<`_=d=cQ3b6Y2Tc8^^`jE/I4k_Ti?/F/0^b
?8`oRHAVZ4ND?8BGBQ9m3?@n8Y?::DQ]Nne9YLbJRDNm2UXc0=?]^/42`<gL]`B]
<7`E<d8UI=j5`[F4;A13Z@XM>=J/=4j>X5?C]AE?To>>^mLAC9RIDTelWDCAdO_D
nk?YgB3h2DcLR=CBe61CKB4Do0L2E7DT/bdFR`@K;mEc9M:h`c?kV^QMd=<5N3U/
O2Q5e]>ER8W55m^eGXi/I@DEX8klGgOO^09CI:TDPXSkd^X^6?AQn:?WR=XREJKd
bN6>>hmK@DhlfL[do0]A:a0mb=i;L_;V/MCUIMBf`V?N3@aNd<ga>SYa6jJnIm4e
MA?PC6FDWDnTiVcVJ99JFd`ebGhKYhaJ:e5XHOZkcHg857a_6dNZOT0@K7T6T6m/
lZ6Z[Z]KFZd/]kU;RG;E6Gc2^FT;EF;UD1NmkCG3]Qj13]l>_`Xm/>aSGhm/H@0a
af@I]P]QD3X<ljK8W<JQRlV_[CN^ak6[?DHX>@Dj<G5lI/O70`ULYZ]D4Zg1L>g[
9UNFEBhMdeQ2nJKo_nLa3WLA90PdCb=8bEhg/GEMGa`8NOeJ=keN>1HAhnVj[R5=
_GJW>hcd=4M9KF?i1KGCNGa_VPaT2WfkeoHV:g5`ICU>hM];Tbk]M?NGOjOP4=AD
iH=BeSD6TDNDXX2UUU><ODleil3`8>[:2EeoQN2@/dA@T@aRAU:Z?D=1mA1Cf?@S
Q8C`aSS87^<I0j]m12VcSNFljhI[lR_;RE:kELJFeF9VPZh5>A6dR8b:0JU@cG7_
XnR>]1JI=2W^YK70meL47X?ZLHG4[AJTD;d;``Q4l@DM6KFHCe/1?V@@TS<34XT3
77A>OB1R=Cjga/M;=2mmbgff6jQ^bYLKmBgZZGC5a@HHW^;CNhT64C7G`8hPm0U@
/S0U?UFe521[jXD]S:NMm:cJR/7OS<]1P:K6P_8Ld]Mch>E=^QT?6MT49V54^OA6
n=TK8b/M9Xfd[K9B/@dI=7T/aH47<]QfFUcbT5;=K9lajE;GV_BTS4XUK:fXgP8V
U9=iccg64j8fG4[/;e`U?MG6n:;cgf8m<g1f8jf2O<KQ045lP:F_:DUL]ZL<4Heb
aMBnRB0:IJbiOQ7LH/oU`^cUhkMo4Q82]ALm<R<SiTJb3A[@YA]3f>cPV[?B=dCR
M;04i;dZJEb^a08c0fBD3XYPHM@22O;6Y6VNDfj=A1?8===9e?V=T;9X=TTPj^M/
WBP=8KC3:P3lbYD2QbgEE@[3=Y1R:XaIc684=OBO>>YTIZ:0HRT11aQ@FNR3;RbL
[d^1RUaNBC2ZlljAbE10X`WBTYRZ@[FiU54GB<dWPJnRV`:HfAVUV36U`@F7W]9L
H`TTfVDE3/;DE<0[3NbXZH3VR2dbPeDK76JE<Ij9_<Z]BF1Y>RB5VZ99T[Y26`dg
MCFFBUVZ:_A8^GMbAh6n:`]KJeYPJRJBb]9?VO=>hBEBZF//X37W@d@XPC`Q1gj_
>d/C>Yl:R;:C7HkHfadR94b_/aEMSdf2;WO/8Z>6A9:C`XRU?=T5Y[<dFPK=U<SB
;e9m><LNmKBKgG[]JglX;S:Cek[gL2gl9;8B>FQ<Jj=O4NeB/O8TQ1JXn2dWKI5S
]nKYB9^G4QIC7YihELLC`TQJWI85=XF0WTe0<CB0AR2TN5ZB@@MFB/]CB5dh[_HV
AM2jVT4T9GLg0<nJn^UTLDmkTA2PCOY`b71/L@Hh<b;26<M4@PV<ZjE:WEUjDFh1
CQ2gZQ[_a>a;?D41:kSF6^ho0aUe0@II_:j_4EIZ[gQ]>10WS</gI_bKHP]AgbJE
kb=P]A5Nf_F<1c2^:UV;LfkFdek7AW[m11`8adUKj9fZ<T<j2oNAm4B@[A0lX928
Q118AW6d1F6TI6a1VS=>9/CJEV/f<X;n9P/aaPFV[aGLG2L98HbIEc]9IlVFNIkV
NdJfSb5P:IfT8T4l^5PnI=SIlB@eZacHbD<oDGSK0lfFGT2dL;2^4/aghViHcENF
/m9M_P@ZU0:`0FSNanC9Zn`k^R24KacX?9=OJe<BP=4O7Q</7K0YE<fj3XJY:</T
22bfTXS@:?=`7<Q2=>ZSZ33DZE]HWkPfk??4FhJ/U09XBPU92<F5]J@kQo1IU_Cn
O`VTPA:l[H?B2ai8?;JEZ6<HY_P>nF_J7iH/gH>0N1j1Q7_MR2d8SGgTSBEaQCk]
]YX2^0ki@MBFYX[NC84KJH4JN8l=PHRfTGJBnKfFd=6D[YGKP2_ZQUUMTF>AWWeC
nV?Zce[KHnOJ6JCWMP?N7L7bl3e`U0WFQmOA_PNb0Q6jWhdN/>^6VMY;`DFjQdnT
_IN<4^G`m8^gFWlR/Zlk0f[KRY9V_NMoU4RMU]NZDD0bLHE6CfW_TC6eUK;N/UE9
2MJdnOBYkm]eT16o;/TOPF0knco4kQPKaG=<4=jXZ93Jg=c8cKad@3eoSejf1Lbi
EcS:_S`]G@=BCBTmJ;kI3k_2Bj7KiDghb9_o_C4]R@Te3AGXY?R5BfFml`JY2:Zd
RhURSaWn<bhG9QLDi0P`^9_1^Ra;[Y:j5l0i/jGaVCT`?HaQd;Ak09d9QhbKk9AJ
QH5;43O;PJE:MKO;A:4JU:AaIA[cZ:1[Rh88UKSh/JC:>XQ44bc25BNCo:lSL:TH
X8T;[8hkJd0YD4147`KcLh]<78BQPCabURUKGa8;NZN0V[2S;ClQlZ=R8c?SH[H4
3C:3mZDaB3P24SfX20X;3;EcJPg:6Bd8PCe2^hQ3ZH<C4Snk<QeJk3Zl0BaD5JOR
CB[kg@JZ5hLd?R5kicbX^8TR9S]d@RYT1l92YY8HWBfcGkI;kkBC0[^BgCPH6<6<
:9/ICCW2^So8X:kJGHGY]:O`UKa;4add=HUK<bJKYcoKKaijcPPSSYU`FSQJO`[D
7d?PLb5h[D8;EQ<kMXRiiC^[jZMRii:3P`nD9b>4A[:C8c`lIJg3kX6;:DAk9lJA
KB/JI/I:ZXWS3AL2Occ0SJF2c45ie:E<=[aCnl=8PD>Jcd>dFQmkiI/lG2KP[1o1
Y1PBg[Vj8k][W:ojULK2^kObB3_3B:_J_;kcXB_38<;;fCFb=RZhWTM[`Kc:Fi:i
];nR9cDCndM3oBZlDSmHYlF;`C`C3PZT?L?Y/h[[TYbI_X?KC;/D@Zh<H[9@c9ZI
m4;T438S8aATLRTKgce3E4P[eCkblPWB[bWjO/N0fPTTSii/Y[mKD<6ShIaPSLRD
CkY:Oi0L7i/j6Q<a_jT18jULIcVdVK@/@h9DUhV4J;U;W[IdXJ<jTLLC:<33]TPL
Q3@D^RBLXc;dD<[A]4N0^<LTSKDdSShJldZ1mASaY4CincJl>`fTTlG<U<NYUZnI
_H@K04HL:[]:PhcXX4;;V/BSH4BdWk<KXTLK;DJ/BKZ:lPZ8ZD_XI/7cc0/3IkY/
Z<dj=/`L=cKXQB[cT273iAClb<7dUHBT8JSCcj93l9:3?4HSBAkSC9aT@Jgc^DOR
VSR:Rd^C4bM/Bb^kL;JLD0ZdGDlD_kTLNPH[]JBJ]3D[DldcDLELmRRcE:71ETa[
nDFlU40kF1oBm;5hDS>RQ2E:iC^CIJhHmPW2E;]C2;lJfl6CJ3<Sa<JLAd;b2@X3
TDN4Fdl[c0[CAM1dLkh4n4EK6[AB^IS/LS>9ECJdhD[Q0<hM5PJY^LGXhQg0]CdI
dCFd/SH3/JAcR043X572oTcjNdSLQlmYeKGdJLPh6:1iTDelN3nB6Cc45de;o:NK
n[mTk58hWcccZP>M3/_8Q@/dTC;e:1n55L7h[/n/h59aKhTTO</[^]0`9DcJM9B1
T48bnZ9djl^M24k2=jM/d/dS8Sg3;b>Z>]38K`km8/k<62jZe4a52X8</4=@V358
S:mXBRkJGD>dGLnP]/b8fRPRC8/`SM7iA<J<3R8Hlc>cjaJM>hlcfRbL_4Z4^=?e
FDB[8<:4ORC9];nDf0fP^<AUAJZ/;ca=H]AZdZLAKCb@AH/3QZKLbJ/8]X[YeXZT
NM5i;=GPZ`]Di;GRZ3eJa]7;dBUIT3>/^dJdBlX4`T^DAd:/:XX;c5<Xfi]8/b=S
=X]Y5E?:=]HdS;;EOQe<SMB<KYDLlLf=J8/4[HKdG5QMQT7aL8^8SJYkP/CDVMRK
Ne824LV[jm>Xm51TJ5O/mDg4J[eeKbIk^H78P820=j1<f2@9_@AV`J4@^28^72V8
2291R:7>;7^<6:=8>>4Z?A@<A8@B2AB:@BB9aBBbVCbd<5:H;fIC>I/FK=FL@J3C
QZ[BOCbM@9cD>RDFQd5_CIRc>O;B<7^=E6Y3VZ46[EB/EV/EHPeR?4Z]eJ8F<DS6
cEbH5:afV^ERcFnhB:ab/@12kGNkgBgfB8FnfeBmgfcGbbCb6<fNHVLDf^BJEB^D
GBBbj@G>9IK9IRGGNWc=FImBJ6UDRRCbI`R3DCDdKFC^LDZVCk@ZBWdm:KN^GlLk
[LeN]FjhL6cERdfnXl3QG6FHk:G:8iO136ah4Di;8Hl@F>UH[^=FWmK9IF8I7VN2
GNG[IB:6[fDgIk?AdR1CbVd2PdAgOWFJFL`:U:DYZV/nESKTY0J>46]j`R2jJc>L
5;NR2hC2951JgRe32^@JjCa^`^Rl;/Ud=/`/K8_8m:FX@`b29hkkX>Zjj@aA6C9X
Y4ZF/a1kd96UZ_`70Y>B2VZK=FLb3?RWCmbDeZP_nVcHUY0LPbVCRX`D`3QZdjT]
^2`Ro;4R2^@`;D[Qc3Jj90^lC@liLfai6Sc>`c4EXBPBW/dmDMc@U[6XT[lm_Dk2
G:a0]2]^Y/T=DnA_9i9M7:>X:4=PVQN_Nd<2XLAJ_a6j:jBl`Jn?560@:`/LaZoD
[2>O=lda2RTeEO42A;bUTLaO7B9>XY3Ie/l<f>cEK8Zi0J;3W0[SZY>;:;NfZ<@:
fJI]8XeYDLo9gDEAKi]O9jI@4clZB312^C;3CXbi3Rb95LR[:n[RY8eE4/GA@3VD
34mP;8gma[>gcM`RZbUEiE@D_>l4H[XZ=<c7:Kf3E6jCgZUZXZOC<Y_ZoY_BI9=Y
Z5Am/>i9bUfjfE;BVYmf]kPJb;Nb;<>X[jXZh[mhZbgCTe05=H9>]l7^YW^K^AOH
HZaH^FCbUMNjCEKO[B[68g/^jgf;RR7@6Y[DHk8ckdOJcl?c8Z58IT<X=W:L2ZO1
NQ>5Cd>hQG</QcQ/a[CVJ?>AWeKXU408KH/Ukh;6X<Ca`BnKkNS[JEDK2Ko@EHPQ
EX8;FV69Xa:L1^fW5/Le9F^bDec^dU9c?=1Bd2f:gCTGEW5kmD^2dgN<ECk]ND7h
OaO6g]YkVL7IDekM=mNl=_?0^IFRj[]Yj@C535Rdc0:OJaSN]HcSbQlkjO?cYT=9
UkKd1UITb<@S2ME`O`7FZiH_I2EVcPm[lZFN=gd:N7PdgD0hb=D14WhGaQkiWL8O
EL1QZ9IWGTJDbP=8HaF`=IJl]MSSe7/:;?/DU9cXFA2mJ/Ci7i]bW]=:doa_c^5>
0YIP^l]RH2b/h>XZ9Nb_fN8_O92/^3ZFgU`;6YU1226R/h?0oAO3_U0>@AfZI31C
g;Y>FXaZ9RT89;I>i1Q/K8d0YBB2YU:cL3Z?oO/F<[2Hda]1:/Um1QdHP>@O>Rm]
A`FI>cGlQ9;YO5V4I8dK<Z<@8AA3EncYDC3CV6j:RPDdi1F]_CF`TE9L3BM7LLnm
/VccdX^T?NYUQ/<gg/4J>_R4hHX@A_Q;9M@3>XAAT=n_XgYb9BUE:]613<:Vn^fS
<oenIO7`U@GPAiV9ET4A96:Xf2<RH7@AT@D>ARfS4aBFhmeilEQ>Af8dQ5E<V4JA
YG`gA^/HR^[bQ5=AU[0iD9JKQ=1NDhPLZhN0bXR1F446c@:EmaSQ9Y5AOG:`;BaI
14bDC4iK3ehVb6PRm55C68<:D9XP=B`Y2VaQQ2oe7ITBg/X22c@gD@7FVKP72:J;
Z9EAU1C=n=jocNcVQ4Q^19WhM@k@BPia3oB__[:NaHWbdIPWkV6XiK<QbP`ECY@>
:U1544f>:_b>0<H0@4N@B]<[LbHD@=nh=QjPC`7DJ55m31GR?CTY1=T[<NGjTBC6
DeC9CgH9mGa35O1Md8UABW0^2Ze95ClWdij@[fF<6VWa<U:5?co`YZ76NC98Rd_[
@RkM6T/9Y?bEnQ]MH@Jn=eIVKaO[hjE6B:h`eQ93SIU?//AY<MI4cY^NF5Y15NX:
^KWmJMj]LGYTlTM1TY]1;Ff`BM=RX[bB:5HWTQRX;@i]>5H>o542XVK`eXfcB`:6
;7A];1LDSc?`HZIX>U>MT7Xj>aUHbZXZ5d<EZB2aLJZfJg42N_N=[laR4>SO15LQ
bQQ:>F<oJmBSg:m?LC:Q^a3_bAFQUB3U]JKWO?6XXRnDeF[811F@3VJck:/`ZRlE
@fJ25WVQBVPQf41`aBFANGL/IDD4D`5XTRVY[:hDd?lCJPf4TPZHDcR/AMJj2:D_
QRmc1?BO5SSF36omAel?[VSQPT/G5F4]T^L:iB2k6eGGHAh]>Al5F<9QLh]=V`a5
LDb/DZ>2ILQ1W=7`RQGl8W_9WJDXf8g[Y7VC@O5@ReVQkbnC<fMN/GTbPEGLoi/d
beVEW0<UdVfQCXH6I5_I;UcW2`6P_@aFZ8BMXPbYLiF2W[5bGD^G5W4<=h11MX;B
6ikhb9^aQJdRW?W^=2YVmT?2;HMB31]0=NJHaDTQS4Q21Lk^[TcQ^51<<O=60aWQ
nk`kQi5GUPTgFO<Tg;`@C2d>S8hUT[5;MVMaijTFedF>Uj@RIHD8eMk6M=I6b?9m
Q>U22<ZjW8i]FIYXMEB@cOZoEjD2U=XXf[eWEP2HF2`H`>Ld0W7YbSUD:oMEPUCd
@FEo1>T@YHl:Y7SOZj66j@dU/lV47bH4P:S8<QIB=]gPP]1U:MHLeTHF:lg5V43K
g^5I^cU4TG_j/@3^mlA_HF^hdaDaNF?WN8a@gAX7=46QBM=eDX:B4N3M29QHYW3:
[;R;j1LYZI62gh^5jb5co6RM5:@7IYQYk2Xl=8^AQ13E15g/G09a0^<MEMWNobOM
Q<bX/dZNNX]3TmWCA/9fnJJ>2nMabGg[H4KmodKkKHk9FTCT:OGA6@YZDgeB/MQV
_=]lDW7M<Ee<X9<i:MJ8eb7[XNm`g/B1ZW]>XQBMY5HK?dU1mEK2:]7ZN6d=5khC
Lh>RcTIdUDYKhcBPF^J0ig/KV_UO>m?[Q;3?3eVKZYS:R/GL=9OBVQPFb631QjfZ
;:Eh[cNh[eZ6k7jRBF9aB85nmmkddTYbN/`/V?XeA:EcBQNW1_AMS2BKn@^m^D[6
JL0IY`CB?nPYG3oBhj/`e68LfJMP34jJ3f^T=2TPnPAHlTRF[TjZ9/K:_DaFlng>
?O0BX>58e8=^^RB30W0^@:OJoDCNP:1RlF2dl:EFEjG^N<_bS6/V;Xh0h<oV:T/j
UF^@b@`8=jbHo@>Rki1P]hKR?I0/b/mbjf?H:NWF=4TL[X^nYV]<7]2TLnX2<B@0
9m0/k4G10]00k2llD:k6UZ@FmTbj8lj0Q/OlKJMg3D9H90XNl6ojI`h<m^me3P`/
fQ1ZXPXT:]3^9PNLI8nn^X3fC6Ob34lBB2WVNH8/_HgF5XnVNTPNF/c:a:PbeD]T
6:P`PjllK<i:aB=^aI50aG4h_JB2`>G:VTH@GN=n_jCNmHcbKhANe/^/S0HJo^^H
2Dh<l2h2^g2>fLXPS8:hdf:K2i32k8=^l/8L:S4Jo6nK0dYTPT_8PP::]I0TI;4m
596<D;6HPFcLD^=^i8llhl3g74=]44R`8e373HE8:Xj0/o5J]nb2Z<gfoLd^;<oh
fP]b2d=dn0:/bc5c5jThl8m/K[10nG52aK0naA0^_>DXkA6JX>:DV0PT?Y40k72n
4Y0i8`nc6a6I2^5Y5756cFn`4G96cF]feh9:`/a^mHVXN<en=jj17bF;3f=bb<hN
b^gZUHkd@Fjn^[0?0?8c8c8LNlm4^W6`c29`NS6X?nTP_@nl@;71514HnLcY:=57
9j3fld3VdN8]5e5>Ag1^[jKJoJR89HTg0?7b`;7a0::W7HUJZddBkU0:h<]`GNoi
2@jbhn8`H[8B5990G11110`m4^I2I0dj6[8lX?8?<G2l4Y4K4>4hg32<34cG1[1D
E694/LkY;6b4ejN>0P;C2Dg;;g06FK1A9Y1ZC7093Nl4]33jbF/La/=jIY1ZnBn^
8i8_8]4k33:VX@Phm8_N9ljTV<:E78FJLY:blfcFi5<K:_=63fjb?HM>f;<^@[;M
7e22L0/8LP;^:`jhi2e=:f^[;`HJJ4^DWY>R3D[hd>f>2T/M;35>>X>8hNo>]d9Q
:18g690_0lBY6I:;6/4hcVVD5XbbR<fR3VWY7c9R8dbU?C9nle5<iPfG;2Rfm[26
3E:dbb[<jG3hhS==@Q;I5dohXPC6bU497]1e9g875`f@g^QX;6oTUI7o0==^46Y>
470c0dD/o4D/Nm6F=2@@ee5C5V`/`ZhNke;G1]4E?5>B3VMQ;@8e097c=]BI?3A3
20XhLT2U51<VTg;`XDSOAW>@Rc24h73PIXQV@g?BWY3a5i<Y=98A@2kB=V^PgA>2
iEAhm6TTo1>05HhM3R1RmM3=?]?LSL/M=O9`j31`2eBNlg>A1@MQ=[81:h3Vo[DH
VO@j2Taj;=7:4YCF@aAQD7D^35?:Gk9/h?CM>VGoE>:o@ZIGA@2E:P8hj0/LjC0]
?hngAeEe;o0_0c74J59X^Doi=Lh<dR@S?OATdSBVmcD5AH2U@Jla83DDZe7B;997
B=D0W8MPV`[C;I?RZ4kLDj>/:hnZB3CnbS2BFM=g2cGFne6_Eg;k2e6bg046@:;B
gYF;3e1e3XXQDBfPkeA73]?Y9c3[5fS1DSCLOGEJ>33@cG7bmg?WD:2EDBVa;[FU
<V^7C<hQEAH>8dBQ8_718_<G:S?cGLBYGU=e8daGC9He?H`Bo`KZ`CEYE53cH5F;
1[AWCIH8`6;d;=A_HO1dXm1WIf:eBAH8[h^6;N8200d1000300000@1V00010@03
00000@2D00010P0300000`00<KH10`0300000@0500011P0300000@0200014@04
00000P00<K`15@0300000@0300015P0400000@0006/15`0400000P00<L@16P05
00000@00<L`16`0500000@00<M@1700300000@010001:00300000@0200000000
00P02008000020008`<002;k000>4@0:o80002L@00[lP0009a001`2^07X01`01
00400`0300@0000=00`0401V06H0^P2j02P0A00F455eJF=[E6U]IJXPHFiT8640
02P0BP095eA9ATHP:4aJEbTPI6ESKfe`LVEcLfmb02P0D0007f5bIB1^IFETIF@P
M6lPLfEU87AXJG<PL6USM7EbIBh0003o\>"], "Picture",
 ImageSize->{102, 148},
 ImageMargins->{{0, 0}, {0, 0}},
 ImageRegion->{{0, 1}, {0, 1}}],

Cell[TextData[{
 ButtonBox["Marshall Rosenbluth",
  BaseStyle->"Hyperlink",
  ButtonData:>{
    URL["http://www.aps.org/praw/nicholso/00winner.cfm"], None}],
 " (1927-2003)"
}], "Text"],

Cell[TextData[{
 "For our atmosphere model, we wish to simulate the random thermal agitation \
of a large collection particles interacting only with gravity (and the \
surface of the Earth.) Basically, each particle makes random steps upward or \
downward. However, these probabilities will no longer be \"unbiased\" (",
 StyleBox["i.e.",
  FontSlant->"Italic"],
 " 50-50 up or down.) Instead, for each potential step, we check the energy \
change involved and reject those steps that are thermally unlikely. Our \
simulation will start at an initial state and continue through a series of \
time steps towards some sort of \"thermal equilibrium\" or a \"relaxed \
state\". \nHow is this done?"
}], "Text"],

Cell["\<\
We will follow a set of specific rules.
Step 1) Initialize the positions of the particles in the atmosphere.
Step 2) Select a random particle. Make a potential random step (up or down).
Step 3) Calculate the energy change to the entire system if the particle \
makes the step.
Step 4a) If the energy change is negative, accept the step.
Step 4b) If the energy change is positive, pick a random number between 0 and \
1. If the random number is less than Exp[ - (energy change)/kT], then accept \
the step.
Step 4c) Otherwise, reject the step and leave the position unchanged.
Step 5) Go to Step 2 and repeat (for a very long time....)\
\>", "Text"]
}, Closed]]
}, Closed]],

Cell[CellGroupData[{

Cell["The Potential Energy of an \"Ideal\" Atmosphere", "Section"],

Cell["\<\
What is the potential energy of our ideal atmosphere? For simplicity, we will \
neglect all interactions between particles. Essentially, the particles feel \
only a gravitational force, and their dynamics are independent of other \
particles.\
\>", "Text"],

Cell[TextData[{
 "The total potential energy for ",
 StyleBox["N",
  FontSlant->"Italic"],
 " particles will be\n\t",
 Cell[BoxData[
  FormBox[
   SubscriptBox["U", "tot"], TraditionalForm]]],
 " = ",
 Cell[BoxData[
  FormBox[
   RowBox[{
    UnderoverscriptBox["\[Sum]", 
     RowBox[{"i", "=", "1"}], "N"], 
    RowBox[{"\[Phi]", "(", 
     SubscriptBox["x", "i"], ")"}]}], TraditionalForm]]],
 "\nwhere \[Phi](",
 Cell[BoxData[
  FormBox[
   SubscriptBox["x", "i"], TraditionalForm]]],
 ") = ",
 StyleBox["m g ",
  FontSlant->"Italic"],
 Cell[BoxData[
  FormBox[
   SubscriptBox["x", "i"], TraditionalForm]]],
 ", is the gravitational potential. Particles are at high potential energy \
when they are high in the sky and lowering potential energy when they are \
close to the ground. (The ground will be at ",
 StyleBox["x",
  FontSlant->"Italic"],
 " = 0.) In order to simulate the ground, we'll take \[Phi] = ",
 Cell[BoxData[
  FormBox[
   SuperscriptBox["10", "5"], TraditionalForm]]],
 " (a huge number) whenever ",
 StyleBox["x",
  FontSlant->"Italic"],
 " < 0."
}], "Text"],

Cell[TextData[{
 "Computation of the change in the atmosphere's energy if a particle takes a \
random step is easy to specify. The change in energy is \[CapitalDelta]",
 StyleBox["U",
  FontSlant->"Italic"],
 "= ",
 StyleBox["m g ",
  FontSlant->"Italic"],
 "\[CapitalDelta]",
 StyleBox["x.",
  FontSlant->"Italic"],
 " If \[CapitalDelta]",
 StyleBox["U",
  FontSlant->"Italic"],
 " > 0, the upward step is accepted if a random number between 0 and 1 is \
less than ",
 StyleBox["Exp[-\[CapitalDelta]U / T]",
  FontWeight->"Bold"],
 "."
}], "Text"],

Cell[TextData[{
 "We characterize thermal equilibrium with a temperature. The energy (per \
degree of freedom, one in our case) of an average particle is ",
 Cell[BoxData[
  FormBox[
   RowBox[{
    SubscriptBox["k", "B"], "T"}], TraditionalForm]]],
 ", where ",
 Cell[BoxData[
  FormBox[
   SubscriptBox["k", "B"], TraditionalForm]]],
 "= 1.3807 \[Cross] ",
 Cell[BoxData[
  FormBox[
   SuperscriptBox["10", 
    RowBox[{"-", "23"}]], TraditionalForm]]],
 " J/K is Boltzman's constant.We will take ",
 Cell[BoxData[
  FormBox[
   SubscriptBox["k", "B"], TraditionalForm]]],
 "=1, so temperature is in energy units. With ",
 StyleBox["m g",
  FontSlant->"Italic"],
 "= 1, the potential energy at ",
 StyleBox["x",
  FontSlant->"Italic"],
 " = 1 is equal to a thermal energy at ",
 StyleBox["T",
  FontSlant->"Italic"],
 " = 1."
}], "Text"]
}, Closed]],

Cell[CellGroupData[{

Cell["A Gravitationally-Biased Random Walker", "Section"],

Cell[TextData[{
 "We now have everything needed to define a \"gravitationally-biased\" random \
walker used to simulate particle confinement within a collisionless \
atmosphere. For simplicity, we set ",
 StyleBox["m",
  FontSlant->"Italic"],
 " = g",
 StyleBox[" = ",
  FontSlant->"Italic"],
 "1."
}], "Text"],

Cell[BoxData[
 RowBox[{
  RowBox[{"s", "[", "]"}], ":=", 
  RowBox[{"\[CapitalDelta]x", " ", 
   RowBox[{"RandomReal", "[", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"-", "1"}], ",", "1"}], "}"}], "]"}]}]}]], "Input",
 CellLabel->"In[1]:="],

Cell[BoxData[
 RowBox[{
  RowBox[{"step", "[", "x_", "]"}], ":=", 
  RowBox[{"Module", "[", 
   RowBox[{
    RowBox[{"{", 
     RowBox[{"trialStep", ",", "\[CapitalDelta]U"}], "}"}], ",", 
    RowBox[{
     RowBox[{"trialStep", "=", 
      RowBox[{"s", "[", "]"}]}], ";", 
     RowBox[{"If", "[", 
      RowBox[{
       RowBox[{
        RowBox[{"x", "+", "trialStep"}], "<", "0"}], ",", 
       RowBox[{
        RowBox[{"reject", "+=", "1"}], ";", 
        RowBox[{"Return", "[", "x", "]"}]}]}], "]"}], ";", 
     RowBox[{"\[CapitalDelta]U", "=", "trialStep"}], ";", 
     RowBox[{"If", "[", 
      RowBox[{
       RowBox[{"\[CapitalDelta]U", "<", "0"}], ",", 
       RowBox[{"x", "+", "trialStep"}], ",", 
       RowBox[{"If", "[", 
        RowBox[{
         RowBox[{
          RowBox[{"RandomReal", "[", "]"}], "<", 
          SuperscriptBox["\[ExponentialE]", 
           RowBox[{"-", 
            FractionBox["\[CapitalDelta]U", "T"]}]]}], ",", 
         RowBox[{"x", "+", "trialStep"}], ",", 
         RowBox[{
          RowBox[{"reject", "+=", "1"}], ";", "x"}]}], "]"}]}], "]"}]}]}], 
   "]"}]}]], "Input",
 CellLabel->"In[2]:="],

Cell[BoxData[{
 RowBox[{
  RowBox[{
   RowBox[{"T", " ", "=", " ", "0.1"}], ";"}], " "}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"\[CapitalDelta]x", " ", "=", " ", "0.1"}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"reject", " ", "=", " ", "0"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{"NestList", "[", 
  RowBox[{"step", ",", " ", "1.0", ",", " ", "10"}], "]"}]}], "Input",
 CellLabel->"In[3]:="],

Cell[BoxData["reject"], "Input",
 CellLabel->"In[7]:="],

Cell[BoxData[{
 RowBox[{
  RowBox[{"reject", " ", "=", " ", "0"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"(", 
   RowBox[{
    RowBox[{"w1", " ", "=", " ", 
     RowBox[{"NestList", "[", 
      RowBox[{"step", ",", "1.0", ",", "50000"}], "]"}]}], ";", 
    RowBox[{"Print", "[", 
     RowBox[{"\"\<Rejection Rate = \>\"", ",", " ", 
      RowBox[{
       RowBox[{"reject", "/", "50000"}], " ", "100.0", " ", "PerCent"}]}], 
     "]"}], ";"}], ")"}], " ", "//", " ", "Timing"}]}], "Input",
 CellLabel->"In[8]:="],

Cell[TextData[{
 "We generally want to have the reject rate be between 30 and 40%. As \
\[CapitalDelta]",
 StyleBox["x",
  FontSlant->"Italic"],
 " increases, then the rejection rate increases. To decrease the rejection \
rate, the value for \[CapitalDelta]",
 StyleBox["x",
  FontSlant->"Italic"],
 " must be decreased. \[CapitalDelta]",
 StyleBox["x",
  FontSlant->"Italic"],
 " can be larger for larger ",
 StyleBox["T",
  FontSlant->"Italic"],
 "."
}], "Text"],

Cell[BoxData[
 RowBox[{"ListPlot", "[", 
  RowBox[{
   RowBox[{"w1", "\[LeftDoubleBracket]", 
    RowBox[{"Range", "[", 
     RowBox[{"1", ",", "50000", ",", "100"}], "]"}], 
    "\[RightDoubleBracket]"}], ",", 
   RowBox[{"AxesLabel", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"\"\<n/100\>\"", ",", "\"\<x\>\""}], "}"}]}]}], "]"}]], "Input",
 CellLabel->"In[10]:="],

Cell[BoxData[
 RowBox[{
  RowBox[{"walkAtmosphere", "[", "n_", "]"}], " ", ":=", " ", 
  RowBox[{"NestList", "[", 
   RowBox[{"step", ",", " ", "1.0", ",", " ", "n"}], "]"}], " ", 
  RowBox[{"(*", " ", 
   StyleBox[
    RowBox[{
     RowBox[{"initial", " ", "height", " ", "is", " ", "x"}], " ", "=", " ", 
     "1.0"}],
    FontColor->RGBColor[0, 1, 0]], " ", "*)"}]}]], "Input",
 CellLabel->"In[11]:="],

Cell["\<\
(Exercise for thought: does the initial condition affect the final thermal \
equilibrium?)\
\>", "Text"],

Cell[CellGroupData[{

Cell["One thousand walkers\[Ellipsis]", "Subsection"],

Cell["We take 250,000 steps\[Ellipsis]", "Text"],

Cell[BoxData[{
 RowBox[{
  RowBox[{"nWalkers", " ", "=", " ", "1000"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"reject", " ", "=", " ", "0"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"walkList", " ", "=", " ", 
   RowBox[{"Table", "[", 
    RowBox[{
     RowBox[{"walkAtmosphere", "[", "250", "]"}], ",", 
     RowBox[{"{", "nWalkers", "}"}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"Print", "[", 
   RowBox[{"\"\<Rejection Rate = \>\"", ",", " ", 
    RowBox[{
     RowBox[{"reject", "/", "250000"}], " ", "100.0", " ", "PerCent"}]}], 
   "]"}], ";"}]}], "Input",
 CellLabel->"In[12]:="],

Cell[BoxData[
 RowBox[{"walkList", " ", "//", " ", "Dimensions", " ", 
  RowBox[{"(*", " ", 
   StyleBox[
    RowBox[{"1000", " ", "walks", " ", "500", " ", "steps", " ", "long"}],
    FontColor->RGBColor[0, 1, 0]], " ", "*)"}]}]], "Input",
 CellLabel->"In[16]:="],

Cell[BoxData[
 RowBox[{"ListAnimate", "[", 
  RowBox[{"Table", "[", 
   RowBox[{
    RowBox[{"ListPlot", "[", 
     RowBox[{
      RowBox[{"walkList", "\[LeftDoubleBracket]", 
       RowBox[{"All", ",", "n"}], "\[RightDoubleBracket]"}], ",", " ", 
      RowBox[{"PlotRange", " ", "\[Rule]", " ", 
       RowBox[{"{", 
        RowBox[{"0", ",", " ", "1.0"}], "}"}]}], ",", " ", 
      RowBox[{"AxesLabel", " ", "\[Rule]", " ", 
       RowBox[{"{", 
        RowBox[{"\"\<walker\>\"", ",", " ", 
         RowBox[{"\"\<x, n= \>\"", " ", "<>", " ", 
          RowBox[{"ToString", "[", 
           RowBox[{"n", "-", "1"}], "]"}]}]}], "}"}]}]}], "]"}], ",", " ", 
    RowBox[{"{", 
     RowBox[{"n", ",", " ", "1", ",", " ", "251", ",", " ", "10"}], "}"}]}], 
   "]"}], "]"}]], "Input",
 CellChangeTimes->{{3.445885234894099*^9, 3.445885246822703*^9}, {
  3.445885334340488*^9, 3.4458853409485283`*^9}},
 CellLabel->"In[39]:="],

Cell[CellGroupData[{

Cell["Statistical Analysis", "Subsubsection"],

Cell[TextData[{
 StyleBox["Mathematica",
  FontSlant->"Italic"],
 " has useful tools for statistical analysis in the \"Add-On\" packages."
}], "Text"],

Cell[BoxData[
 RowBox[{
  RowBox[{"(", 
   RowBox[{
    RowBox[{"Needs", "[", "\"\<BarCharts`\>\"", "]"}], ";", 
    RowBox[{"Needs", "[", "\"\<Histograms`\>\"", "]"}], ";", 
    RowBox[{"Needs", "[", "\"\<PieCharts`\>\"", "]"}]}], ")"}], 
  ";"}]], "Input",
 CellChangeTimes->{3.445885358881974*^9},
 CellLabel->"In[40]:="]
}, Closed]],

Cell[CellGroupData[{

Cell["BinCounts", "Subsubsection"],

Cell[BoxData[
 RowBox[{"BinCounts", "[", 
  RowBox[{
   RowBox[{"walkList", "\[LeftDoubleBracket]", 
    RowBox[{"All", ",", "1"}], "\[RightDoubleBracket]"}], ",", " ", 
   RowBox[{"{", 
    RowBox[{"0.0", ",", " ", "1.5", ",", "0.05"}], "}"}]}], "]"}]], "Input",
 CellLabel->"In[20]:="],

Cell[BoxData[
 RowBox[{"BinCounts", "[", 
  RowBox[{
   RowBox[{"walkList", "\[LeftDoubleBracket]", 
    RowBox[{"All", ",", "100"}], "\[RightDoubleBracket]"}], ",", " ", 
   RowBox[{"{", 
    RowBox[{"0.0", ",", " ", "1.5", ",", "0.05"}], "}"}]}], "]"}]], "Input",
 CellLabel->"In[21]:="],

Cell[BoxData[
 RowBox[{
  RowBox[{"grid", " ", "=", " ", 
   RowBox[{"Table", "[", 
    RowBox[{
     RowBox[{"BinCounts", "[", 
      RowBox[{
       RowBox[{"walkList", "\[LeftDoubleBracket]", 
        RowBox[{"All", ",", "n"}], "\[RightDoubleBracket]"}], ",", " ", 
       RowBox[{"{", 
        RowBox[{"0.0", ",", " ", "1.5", ",", "0.05"}], "}"}]}], "]"}], ",", 
     " ", 
     RowBox[{"{", 
      RowBox[{"n", ",", " ", "1", ",", " ", "251"}], "}"}]}], "]"}]}], 
  ";"}]], "Input",
 CellLabel->"In[22]:="],

Cell[BoxData[
 RowBox[{"ListDensityPlot", "[", 
  RowBox[{"grid", ",", 
   RowBox[{"Mesh", "\[Rule]", "False"}], ",", 
   RowBox[{"ColorFunction", "\[Rule]", "Hue"}], ",", 
   RowBox[{"PlotRange", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"1", ",", "100"}], "}"}]}], ",", 
   RowBox[{"FrameLabel", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"\"\<Box(x/\[CapitalDelta]x), \>\"", "<>", 
       RowBox[{"ToString", "[", "nWalkers", "]"}], "<>", "\"\< Walkers\>\""}],
       ",", "\"\<n\>\""}], "}"}]}]}], "]"}]], "Input",
 CellLabel->"In[23]:="],

Cell[BoxData[{
 RowBox[{
  RowBox[{"p1", "=", 
   RowBox[{"ListPlot", "[", 
    RowBox[{
     RowBox[{"grid", "\[LeftDoubleBracket]", 
      RowBox[{"15", ",", "All"}], "\[RightDoubleBracket]"}], ",", 
     RowBox[{"PlotRange", "\[Rule]", "All"}], ",", 
     RowBox[{"PlotStyle", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{
        RowBox[{"RGBColor", "[", 
         RowBox[{"1", ",", "0", ",", "0"}], "]"}], ",", 
        RowBox[{"PointSize", "[", "0.01`", "]"}]}], "}"}]}], ",", 
     RowBox[{"DisplayFunction", "\[Rule]", "Identity"}]}], "]"}]}], 
  ";"}], "\n", 
 RowBox[{
  RowBox[{"p2", "=", 
   RowBox[{"ListPlot", "[", 
    RowBox[{
     RowBox[{"grid", "\[LeftDoubleBracket]", 
      RowBox[{"200", ",", "All"}], "\[RightDoubleBracket]"}], ",", 
     RowBox[{"PlotRange", "\[Rule]", "All"}], ",", 
     RowBox[{"PlotStyle", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{
        RowBox[{"RGBColor", "[", 
         RowBox[{"0", ",", "0", ",", "1"}], "]"}], ",", 
        RowBox[{"PointSize", "[", "0.01`", "]"}]}], "}"}]}], ",", 
     RowBox[{"DisplayFunction", "\[Rule]", "Identity"}]}], "]"}]}], 
  ";"}], "\n", 
 RowBox[{"Show", "[", 
  RowBox[{"p1", ",", "p2", ",", 
   RowBox[{"DisplayFunction", "\[Rule]", "$DisplayFunction"}], ",", 
   RowBox[{"AxesOrigin", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"0", ",", "0"}], "}"}]}], ",", 
   RowBox[{"AxesLabel", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"\"\<Box(x)\>\"", ",", "\"\<Count\>\""}], "}"}]}]}], 
  "]"}]}], "Input",
 CellLabel->"In[24]:="]
}, Closed]],

Cell[CellGroupData[{

Cell["Histogram", "Subsubsection"],

Cell[BoxData[
 RowBox[{"p3", "=", 
  RowBox[{"Histogram", "[", 
   RowBox[{
    RowBox[{"walkList", "\[LeftDoubleBracket]", 
     RowBox[{"All", ",", "100"}], "\[RightDoubleBracket]"}], ",", 
    RowBox[{"HistogramRange", "\[Rule]", 
     RowBox[{"{", 
      RowBox[{"0", ",", "1.5`"}], "}"}]}], ",", 
    RowBox[{"HistogramCategories", "\[Rule]", "45"}], ",", 
    RowBox[{"PlotLabel", "\[Rule]", 
     RowBox[{"\"\<h[j,n], n= \>\"", "<>", 
      RowBox[{"ToString", "[", "100", "]"}]}]}], ",", 
    RowBox[{"AxesLabel", "\[Rule]", 
     RowBox[{"{", 
      RowBox[{"\"\<x\>\"", ",", "\"\<\>\""}], "}"}]}]}], "]"}]}]], "Input",
 CellLabel->"In[27]:="],

Cell[BoxData[
 RowBox[{
  RowBox[{"histograms", " ", "=", " ", 
   RowBox[{"Table", "[", 
    RowBox[{
     RowBox[{
      RowBox[{"p1", " ", "=", " ", 
       RowBox[{"Histogram", "[", 
        RowBox[{
         RowBox[{"walkList", "\[LeftDoubleBracket]", 
          RowBox[{"All", ",", "n"}], "\[RightDoubleBracket]"}], ",", " ", 
         RowBox[{"HistogramRange", " ", "\[Rule]", " ", 
          RowBox[{"{", 
           RowBox[{"0", ",", " ", "1.5"}], "}"}]}], ",", 
         RowBox[{"HistogramCategories", " ", "\[Rule]", " ", "45"}], ",", 
         RowBox[{"PlotLabel", " ", "\[Rule]", " ", 
          RowBox[{"\"\<h[j,n], n= \>\"", " ", "<>", " ", 
           RowBox[{"ToString", "[", "n", "]"}]}]}], ",", " ", 
         RowBox[{"AxesLabel", " ", "\[Rule]", " ", 
          RowBox[{"{", 
           RowBox[{"\"\<x/\[CapitalDelta]x\>\"", ",", "\"\<\>\""}], "}"}]}], 
         ",", " ", 
         RowBox[{"DisplayFunction", " ", "\[Rule]", " ", "Identity"}]}], 
        "]"}]}], ";", "\[IndentingNewLine]", 
      RowBox[{"Show", "[", 
       RowBox[{"p1", ",", " ", 
        RowBox[{"PlotRange", " ", "\[Rule]", " ", 
         RowBox[{"{", 
          RowBox[{"All", ",", 
           RowBox[{"{", 
            RowBox[{"0", ",", " ", "150"}], "}"}]}], "}"}]}], ",", 
        RowBox[{"AxesOrigin", " ", "\[Rule]", " ", 
         RowBox[{"{", 
          RowBox[{"0", ",", "0"}], "}"}]}], ",", " ", 
        RowBox[{
        "DisplayFunction", " ", "\[Rule]", " ", "$DisplayFunction"}]}], 
       "]"}]}], ",", " ", 
     RowBox[{"{", 
      RowBox[{"n", ",", " ", "1", ",", " ", "251", ",", " ", "10"}], "}"}]}], 
    "]"}]}], ";"}]], "Input",
 CellLabel->"In[28]:="],

Cell[BoxData[
 RowBox[{"ListAnimate", "[", "histograms", "]"}]], "Input",
 CellChangeTimes->{{3.445885392348763*^9, 3.445885400035577*^9}},
 CellLabel->"In[41]:="]
}, Closed]],

Cell[CellGroupData[{

Cell["Mean and Standard Deviation", "Subsubsection"],

Cell[BoxData[{
 RowBox[{
  RowBox[{"avgX", "=", 
   RowBox[{"Mean", "/@", 
    RowBox[{"Transpose", "[", "walkList", "]"}]}]}], ";"}], "\n", 
 RowBox[{"avgPlot1", "=", 
  RowBox[{"ListPlot", "[", 
   RowBox[{"avgX", ",", 
    RowBox[{"AxesLabel", "\[Rule]", 
     RowBox[{"{", 
      RowBox[{"\"\<n\>\"", ",", "\"\<<x>\>\""}], "}"}]}]}], "]"}]}]}], "Input",\

 CellLabel->"In[29]:="],

Cell[BoxData[{
 RowBox[{
  RowBox[{"varX", "=", 
   RowBox[{"Variance", "/@", 
    RowBox[{"Transpose", "[", "walkList", "]"}]}]}], ";"}], "\n", 
 RowBox[{"ListPlot", "[", 
  RowBox[{"varX", ",", 
   RowBox[{"AxesLabel", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"\"\<n\>\"", ",", "\"\<Var(x)\>\""}], "}"}]}]}], "]"}]}], "Input",\

 CellLabel->"In[31]:="],

Cell[BoxData[{
 RowBox[{
  RowBox[{"stdDevX", "=", 
   RowBox[{"StandardDeviation", "/@", 
    RowBox[{"Transpose", "[", "walkList", "]"}]}]}], ";"}], "\n", 
 RowBox[{"ListPlot", "[", 
  RowBox[{"stdDevX", ",", 
   RowBox[{"AxesLabel", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"\"\<n\>\"", ",", "\"\<StdDev(x)\>\""}], "}"}]}]}], 
  "]"}]}], "Input",
 CellLabel->"In[33]:="],

Cell[TextData[{
 "The standard deviation indicates the magnitude of the \"fluctuations\" in \
",
 StyleBox["x.",
  FontSlant->"Italic"]
}], "Text"]
}, Closed]],

Cell[CellGroupData[{

Cell["Entropy", "Subsubsection"],

Cell[TextData[{
 "Entropy is a measure of the disorder of the distribution of random walkers. \
For particles with negligible kinetic energy, entropy is defined simply as ",
 StyleBox["S",
  FontSlant->"Italic"],
 "(",
 StyleBox["t",
  FontSlant->"Italic"],
 ") = \[Dash] ",
 Cell[BoxData[
  FormBox[
   RowBox[{"\[Integral]", 
    RowBox[{
     RowBox[{"P", "(", 
      RowBox[{"x", ",", "t"}], ")"}], "  ", 
     RowBox[{"Log", "[", 
      RowBox[{"P", "(", 
       RowBox[{"x", ",", "t"}], ")"}], "]"}], 
     RowBox[{"\[DifferentialD]", "x"}]}]}], TraditionalForm]]],
 ", where ",
 StyleBox["P",
  FontSlant->"Italic"],
 "(",
 StyleBox["x,t",
  FontSlant->"Italic"],
 ") is the probability of a walker being at cell ",
 StyleBox["x",
  FontSlant->"Italic"],
 " at time ",
 StyleBox["t",
  FontSlant->"Italic"],
 ". Initially, all of the walkers are located at ",
 StyleBox["x",
  FontSlant->"Italic"],
 " = 1 with probability of 100%, The entropy is zero. As the probability \
become more uniform in ",
 StyleBox["x",
  FontSlant->"Italic"],
 ", then the entropy will increase. As the particles accelerate to  the \
ground, the entropy ",
 StyleBox["appears",
  FontSlant->"Italic"],
 " to increase and then",
 " decrease. In a ",
 StyleBox["non-equilibrium ",
  FontSlant->"Italic"],
 "transient state of our simulation, our definition is not very useful.",
 " Recall, our estimate of the probability within a gird-box is the ratio, ",
 StyleBox["P",
  FontSlant->"Italic"],
 "(",
 StyleBox["x",
  FontSlant->"Italic"],
 ",",
 StyleBox["t",
  FontSlant->"Italic"],
 ") \[TildeTilde] (# of walkers within grid box ",
 StyleBox["x",
  FontSlant->"Italic"],
 ")/(total number of walkers)."
}], "Text"],

Cell[BoxData[
 RowBox[{"Dimensions", "[", "grid", "]"}]], "Input",
 CellLabel->"In[35]:="],

Cell[BoxData[
 RowBox[{
  RowBox[{"entropyGrid", " ", "=", " ", 
   RowBox[{"Table", "[", 
    RowBox[{
     RowBox[{"Plus", " ", "@@", " ", 
      RowBox[{"(", 
       RowBox[{
        RowBox[{
         RowBox[{"If", "[", 
          RowBox[{
           RowBox[{"#", " ", ">", " ", "1"}], ",", 
           RowBox[{"-", 
            RowBox[{"N", "[", 
             RowBox[{
              RowBox[{"(", 
               RowBox[{"#", "/", "nWalkers"}], ")"}], " ", 
              RowBox[{"Log", "[", 
               RowBox[{"(", 
                RowBox[{"#", "/", "nWalkers"}], ")"}], "]"}]}], "]"}]}], ",", 
           "0"}], "]"}], "&"}], " ", "/@", "  ", 
        RowBox[{"grid", "\[LeftDoubleBracket]", 
         RowBox[{"t", ",", "All"}], "\[RightDoubleBracket]"}]}], ")"}]}], ",",
      " ", 
     RowBox[{"{", 
      RowBox[{"t", ",", " ", "1", ",", " ", "250"}], "}"}]}], "]"}]}], 
  ";"}]], "Input",
 CellLabel->"In[36]:="],

Cell[BoxData[
 RowBox[{"ListPlot", "[", 
  RowBox[{"entropyGrid", ",", 
   RowBox[{"AxesLabel", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"\"\<n\>\"", ",", 
      RowBox[{"\"\<Entropy[n], \>\"", "<>", 
       RowBox[{"ToString", "[", "nWalkers", "]"}], "<>", 
       "\"\< Walkers\>\""}]}], "}"}]}], ",", 
   RowBox[{"Joined", "\[Rule]", "True"}], ",", 
   RowBox[{"PlotRange", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"0", ",", "4"}], "}"}]}]}], "]"}]], "Input",
 CellLabel->"In[37]:="]
}, Closed]]
}, Closed]]
}, Closed]],

Cell[CellGroupData[{

Cell["Potential Energy", "Section"],

Cell[TextData[{
 "The average potential energy of the particles within our model atmosphere \
shows relaxation towards lower energy as thermal equilibrium is approached. \
For our simple model, \[LeftAngleBracket]",
 StyleBox["U",
  FontSlant->"Italic"],
 "\[RightAngleBracket] ",
 StyleBox["=",
  FontVariations->{"CompatibilityType"->0}],
 " ",
 Cell[BoxData[
  FormBox[
   RowBox[{
    UnderoverscriptBox["\[Sum]", 
     RowBox[{"i", "=", "1"}], "nWalkers"], 
    RowBox[{"\[Phi]", "(", 
     SubscriptBox["x", "i"], ")"}]}], TraditionalForm]]],
 " = ",
 StyleBox["m g ",
  FontSlant->"Italic"],
 "\[LeftAngleBracket]",
 StyleBox["x",
  FontSlant->"Italic"],
 "\[RightAngleBracket]."
}], "Text"],

Cell[BoxData[{
 RowBox[{
  RowBox[{"avgU", "=", "avgX"}], ";"}], "\n", 
 RowBox[{"ListPlot", "[", 
  RowBox[{"avgU", ",", 
   RowBox[{"AxesLabel", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"\"\<n\>\"", ",", "\"\<Avg Potential Energy\>\""}], "}"}]}]}], 
  "]"}]}], "Input",
 CellLabel->"In[42]:="],

Cell["\<\
In thermal equilibrium, the average potential energy should equal the \
temperature\[Ellipsis]\
\>", "Text"],

Cell[BoxData[
 RowBox[{"equilU", " ", "=", " ", 
  RowBox[{"Mean", "[", 
   RowBox[{"avgX", "\[LeftDoubleBracket]", 
    RowBox[{"Range", "[", 
     RowBox[{"200", ",", "250"}], "]"}], "\[RightDoubleBracket]"}], 
   "]"}]}]], "Input",
 CellLabel->"In[44]:="],

Cell[TextData[{
 "Using the Boltzman distribution for the probability function, the \
equilibrium value of ",
 StyleBox["x",
  FontSlant->"Italic"],
 " is proportional to ",
 Cell[BoxData[
  FormBox[
   RowBox[{"\[Integral]", 
    RowBox[{"x", " ", 
     RowBox[{"Exp", "[", 
      RowBox[{
       RowBox[{"-", 
        RowBox[{"U", "(", "x", ")"}]}], "/", "T"}], "]"}], 
     RowBox[{"\[DifferentialD]", "x"}]}]}], TraditionalForm]]],
 "."
}], "Text"],

Cell[BoxData[
 RowBox[{
  RowBox[{"Integrate", "[", 
   RowBox[{
    RowBox[{"x", " ", 
     RowBox[{"Exp", "[", 
      RowBox[{
       RowBox[{"-", "x"}], "/", "T"}], "]"}]}], ",", " ", 
    RowBox[{"{", 
     RowBox[{"x", ",", " ", "0", ",", " ", "Infinity"}], "}"}]}], "]"}], "/", 
  
  RowBox[{"Integrate", "[", 
   RowBox[{
    RowBox[{"Exp", "[", 
     RowBox[{
      RowBox[{"-", "x"}], "/", "T"}], "]"}], ",", " ", 
    RowBox[{"{", 
     RowBox[{"x", ",", " ", "0", ",", " ", "Infinity"}], "}"}]}], 
   "]"}]}]], "Input",
 CellLabel->"In[45]:="],

Cell[TextData[{
 "The mean fluctuation in the potential energy in time is ",
 Cell[BoxData[
  FormBox[
   SqrtBox[
    RowBox[{"\[LeftAngleBracket]", 
     SuperscriptBox[
      RowBox[{"(", 
       RowBox[{
        SubscriptBox["U", "n"], " ", "-", " ", 
        SubscriptBox["U", "avg"]}], ")"}], "2"], "\[RightAngleBracket]"}]], 
   TraditionalForm]]],
 ", or"
}], "Text"],

Cell[BoxData[
 RowBox[{"avg\[Delta]U", " ", "=", " ", 
  RowBox[{"Sqrt", "[", 
   RowBox[{"Mean", "[", 
    RowBox[{"Table", "[", 
     RowBox[{
      RowBox[{
       RowBox[{"(", 
        RowBox[{
         RowBox[{
         "avgU", "\[LeftDoubleBracket]", "n", "\[RightDoubleBracket]"}], "  ",
          "-", " ", "equilU"}], ")"}], "^", "2"}], ",", " ", 
      RowBox[{"{", 
       RowBox[{"n", ",", " ", "200", ",", " ", "250"}], "}"}]}], "]"}], "]"}],
    "]"}]}]], "Input",
 CellLabel->"In[46]:="],

Cell[TextData[{
 "And this gives a fluctuation about 1/40 the size of \[LeftAngleBracket]",
 StyleBox["U",
  FontSlant->"Italic"],
 "\[RightAngleBracket] (about \[PlusMinus] 2.5 %). "
}], "Text"]
}, Closed]],

Cell[CellGroupData[{

Cell["Convergence in Monte Carlo Calculations", "Section"],

Cell[TextData[{
 "Perhaps the more important of the practical questions concerning Monte \
Carlo calculations is one of convergence. How many particles are required to \
estimate accurate averages? The answer to this question can be found in what \
we already know about random walks. The mean square distance from the average \
increases linearly with the number of steps in a random process. If we are \
trying to estimate the ",
 StyleBox["true",
  FontSlant->"Italic"],
 " equilibrium potential energy (",
 Cell[BoxData[
  FormBox[
   SubscriptBox["U", "true"], TraditionalForm]]],
 ")",
 " of a system by averaging over ",
 StyleBox["N",
  FontSlant->"Italic"],
 " steps, then the calculated average would be \[LeftAngleBracket]",
 StyleBox["U",
  FontSlant->"Italic"],
 "\[RightAngleBracket] ~ ( ",
 StyleBox["N ",
  FontSlant->"Italic"],
 Cell[BoxData[
  FormBox[
   SubscriptBox["U", "true"], TraditionalForm]]],
 StyleBox[" + ",
  FontSlant->"Italic"],
 Cell[BoxData[
  FormBox[
   SqrtBox[
    RowBox[{"N", " ", 
     RowBox[{"\[LeftAngleBracket]", 
      SuperscriptBox[
       StyleBox[
        RowBox[{"\[Delta]", 
         StyleBox["U",
          FontSlant->"Italic"]}]], "2"], "\[RightAngleBracket]"}]}]], 
   TraditionalForm]]],
 ")",
 " /",
 StyleBox["N.",
  FontSlant->"Italic"],
 " The term ",
 Cell[BoxData[
  FormBox[
   SqrtBox[
    RowBox[{"N", " ", 
     RowBox[{"\[LeftAngleBracket]", 
      SuperscriptBox[
       StyleBox[
        RowBox[{"\[Delta]", 
         StyleBox["U",
          FontSlant->"Italic"]}]], "2"], "\[RightAngleBracket]"}]}]], 
   TraditionalForm]]],
 " ~ \[CapitalDelta]",
 StyleBox["U",
  FontSlant->"Italic"],
 " is the deviation of the potential energy due to random diffusion. The \
value of the average is \[LeftAngleBracket]",
 StyleBox["U",
  FontSlant->"Italic"],
 "\[RightAngleBracket] ~ ",
 Cell[BoxData[
  FormBox[
   SubscriptBox["U", "true"], TraditionalForm]]],
 StyleBox[" + ",
  FontSlant->"Italic"],
 Cell[BoxData[
  FormBox[
   SqrtBox[
    RowBox[{
     RowBox[{"\[LeftAngleBracket]", 
      SuperscriptBox[
       StyleBox[
        RowBox[{"\[Delta]", 
         StyleBox["U",
          FontSlant->"Italic"]}]], "2"], "\[RightAngleBracket]"}], "/", 
     "N"}]], TraditionalForm]]],
 ". As ",
 StyleBox["N",
  FontSlant->"Italic"],
 " \[Rule] \[Infinity], the difference (",
 StyleBox["i.e.",
  FontSlant->"Italic"],
 " the error) between the Monte Carlo estimate of the error and the ",
 StyleBox["true",
  FontSlant->"Italic"],
 " equilibrium potential energy decreases as 1/",
 Cell[BoxData[
  FormBox[
   SqrtBox["N"], TraditionalForm]]],
 ". "
}], "Text"],

Cell["\<\
For our first calculation, we had 1000 walkers.If we had only 62 walkers (16 \
times fewer), the error in our estimate of \[LeftAngleBracket]U\
\[RightAngleBracket] will be four (4) times larger.\
\>", "Text"],

Cell[BoxData[{
 RowBox[{
  RowBox[{
   RowBox[{"T", " ", "=", " ", "0.1"}], ";"}], " "}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"\[CapitalDelta]x", " ", "=", " ", "0.1"}], ";"}]}], "Input",
 CellLabel->"In[47]:="],

Cell[BoxData[{
 RowBox[{
  RowBox[{"nWalkers", " ", "=", " ", "62"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"reject", " ", "=", " ", "0"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"walkList", " ", "=", " ", 
   RowBox[{"Table", "[", 
    RowBox[{
     RowBox[{"walkAtmosphere", "[", "250", "]"}], ",", 
     RowBox[{"{", "nWalkers", "}"}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"Print", "[", 
   RowBox[{"\"\<Rejection Rate = \>\"", ",", " ", 
    RowBox[{
     RowBox[{"reject", "/", 
      RowBox[{"(", 
       RowBox[{"250", " ", "62"}], ")"}]}], " ", "100.0", " ", "PerCent"}]}], 
   "]"}], ";"}]}], "Input",
 CellLabel->"In[49]:="],

Cell[BoxData[{
 RowBox[{
  RowBox[{"avgX", "=", 
   RowBox[{"Mean", "/@", 
    RowBox[{"Transpose", "[", "walkList", "]"}]}]}], ";"}], "\n", 
 RowBox[{"avgPlot2", "=", 
  RowBox[{"ListPlot", "[", 
   RowBox[{"avgX", ",", 
    RowBox[{"AxesLabel", "\[Rule]", 
     RowBox[{"{", 
      RowBox[{"\"\<n\>\"", ",", "\"\<<x>\>\""}], "}"}]}], ",", 
    RowBox[{"PlotRange", "\[Rule]", "All"}], ",", 
    RowBox[{"AxesOrigin", "\[Rule]", 
     RowBox[{"{", 
      RowBox[{"0", ",", "0"}], "}"}]}], ",", 
    RowBox[{"PlotStyle", "\[Rule]", 
     RowBox[{"{", 
      RowBox[{"RGBColor", "[", 
       RowBox[{"1", ",", "0", ",", "0"}], "]"}], "}"}]}]}], "]"}]}]}], "Input",\

 CellLabel->"In[53]:="],

Cell[BoxData[
 RowBox[{"equilU", " ", "=", " ", 
  RowBox[{"Mean", "[", 
   RowBox[{"avgX", "\[LeftDoubleBracket]", 
    RowBox[{"Range", "[", 
     RowBox[{"200", ",", "250"}], "]"}], "\[RightDoubleBracket]"}], 
   "]"}]}]], "Input",
 CellLabel->"In[55]:="],

Cell[BoxData[
 RowBox[{
  RowBox[{"avgU", " ", "=", " ", "avgX"}], ";"}]], "Input",
 CellLabel->"In[56]:="],

Cell[BoxData[
 RowBox[{"avg\[Delta]U", " ", "=", " ", 
  RowBox[{"Sqrt", "[", 
   RowBox[{"Mean", "[", 
    RowBox[{"Table", "[", 
     RowBox[{
      RowBox[{
       RowBox[{"(", 
        RowBox[{
         RowBox[{
         "avgU", "\[LeftDoubleBracket]", "n", "\[RightDoubleBracket]"}], "  ",
          "-", " ", "equilU"}], ")"}], "^", "2"}], ",", " ", 
      RowBox[{"{", 
       RowBox[{"n", ",", " ", "200", ",", " ", "250"}], "}"}]}], "]"}], "]"}],
    "]"}]}]], "Input",
 CellLabel->"In[57]:="],

Cell[TextData[{
 "This uncertainty (\[PlusMinus] 8.3 %) is about 4 times larger than with \
1000 ideal atmospheric particles, as expected. The reality of Monte Carlo \
calculations is the need for a",
 StyleBox[" large number of realizations",
  FontSlant->"Italic"],
 " with which to find averages!"
}], "Text"],

Cell[BoxData[
 RowBox[{"Show", "[", 
  RowBox[{"avgPlot1", ",", "avgPlot2", ",", 
   RowBox[{"PlotLabel", "\[Rule]", "\"\<62 vs 1000 Walkers\>\""}]}], 
  "]"}]], "Input",
 CellLabel->"In[58]:="]
}, Closed]],

Cell[CellGroupData[{

Cell["Summary", "Section"],

Cell[TextData[{
 "The Rosenbluth-Teller-Metropolis algorithm is a general prescription for \
finding properties of systems at thermal equilibrium. Particles, or \
components of the larger system, are allowed to make random (thermal motions) \
and these motions are tested using a Boltzman probability distribution \
function. If the motion causes the system energy to increase, then the motion \
is rejected if the change in energy is too large. In this way, the \
\"simulated\" random motions of the system eventually relax to a state of \
thermal equilibrium. The properties of the thermal equilibrium can be \
estimated by taking averages of the simulation. The accuracy of the averages \
(as compared to the \"true\" equilibrium properties) improves as 1/",
 Cell[BoxData[
  FormBox[
   SqrtBox["N"], TraditionalForm]]],
 "when the number of particles (and configuration steps) increase. "
}], "Text"]
}, Closed]]
}, Open  ]]
},
WindowToolbars->{},
WindowSize->{746, 777},
WindowMargins->{{432, Automatic}, {Automatic, 23}},
FrontEndVersion->"6.0 for Mac OS X x86 (32-bit) (May 21, 2008)",
StyleDefinitions->"TutorialBook.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[CellGroupData[{
Cell[590, 23, 36, 0, 51, "Title"],
Cell[629, 25, 60, 3, 41, "Subsubtitle"],
Cell[CellGroupData[{
Cell[714, 32, 31, 0, 86, "Section"],
Cell[748, 34, 214, 4, 41, "Text"],
Cell[CellGroupData[{
Cell[987, 42, 41, 0, 34, "Subsection"],
Cell[1031, 44, 762, 15, 114, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[1830, 64, 71, 0, 34, "Subsection"],
Cell[1904, 66, 535, 13, 77, "Text"],
Cell[2442, 81, 895, 13, 149, "Text"],
Cell[3340, 96, 17978, 278, 156, "Picture"],
Cell[21321, 376, 185, 6, 23, "Text"],
Cell[21509, 384, 706, 12, 125, "Text"],
Cell[22218, 398, 658, 12, 251, "Text"]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[22925, 416, 66, 0, 54, "Section"],
Cell[22994, 418, 266, 5, 41, "Text"],
Cell[23263, 425, 1083, 39, 125, "Text"],
Cell[24349, 466, 548, 19, 60, "Text"],
Cell[24900, 487, 839, 30, 60, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[25776, 522, 57, 0, 54, "Section"],
Cell[25836, 524, 310, 10, 41, "Text"],
Cell[26149, 536, 245, 8, 24, "Input"],
Cell[26397, 546, 1136, 33, 130, "Input"],
Cell[27536, 581, 415, 11, 78, "Input"],
Cell[27954, 594, 55, 1, 26, "Input"],
Cell[28012, 597, 524, 14, 81, "Input"],
Cell[28539, 613, 464, 16, 42, "Text"],
Cell[29006, 631, 367, 10, 26, "Input"],
Cell[29376, 643, 404, 11, 46, "Input"],
Cell[29783, 656, 114, 3, 23, "Text"],
Cell[CellGroupData[{
Cell[29922, 663, 53, 0, 34, "Subsection"],
Cell[29978, 665, 48, 0, 24, "Text"],
Cell[30029, 667, 636, 18, 79, "Input"],
Cell[30668, 687, 264, 6, 26, "Input"],
Cell[30935, 695, 920, 22, 80, "Input"],
Cell[CellGroupData[{
Cell[31880, 721, 45, 0, 33, "Subsubsection"],
Cell[31928, 723, 150, 4, 23, "Text"],
Cell[32081, 729, 324, 9, 26, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[32442, 743, 34, 0, 33, "Subsubsection"],
Cell[32479, 745, 287, 7, 26, "Input"],
Cell[32769, 754, 289, 7, 26, "Input"],
Cell[33061, 763, 511, 15, 44, "Input"],
Cell[33575, 780, 558, 14, 79, "Input"],
Cell[34136, 796, 1517, 41, 153, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[35690, 842, 34, 0, 33, "Subsubsection"],
Cell[35727, 844, 652, 16, 63, "Input"],
Cell[36382, 862, 1670, 40, 135, "Input"],
Cell[38055, 904, 163, 3, 26, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[38255, 912, 52, 0, 33, "Subsubsection"],
Cell[38310, 914, 383, 12, 45, "Input"],
Cell[38696, 928, 356, 11, 45, "Input"],
Cell[39055, 941, 375, 11, 45, "Input"],
Cell[39433, 954, 147, 5, 23, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[39617, 964, 32, 0, 33, "Subsubsection"],
Cell[39652, 966, 1703, 60, 138, "Text"],
Cell[41358, 1028, 90, 2, 26, "Input"],
Cell[41451, 1032, 927, 27, 97, "Input"],
Cell[42381, 1061, 490, 13, 80, "Input"]
}, Closed]]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[42932, 1081, 35, 0, 54, "Section"],
Cell[42970, 1083, 698, 24, 47, "Text"],
Cell[43671, 1109, 297, 9, 43, "Input"],
Cell[43971, 1120, 118, 3, 24, "Text"],
Cell[44092, 1125, 258, 7, 26, "Input"],
Cell[44353, 1134, 452, 16, 47, "Text"],
Cell[44808, 1152, 554, 19, 45, "Input"],
Cell[45365, 1173, 375, 13, 35, "Text"],
Cell[45743, 1188, 502, 15, 26, "Input"],
Cell[46248, 1205, 195, 5, 24, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[46480, 1215, 58, 0, 54, "Section"],
Cell[46541, 1217, 2625, 93, 180, "Text"],
Cell[49169, 1312, 219, 4, 42, "Text"],
Cell[49391, 1318, 219, 6, 42, "Input"],
Cell[49613, 1326, 687, 20, 79, "Input"],
Cell[50303, 1348, 689, 20, 80, "Input"],
Cell[50995, 1370, 258, 7, 26, "Input"],
Cell[51256, 1379, 107, 3, 24, "Input"],
Cell[51366, 1384, 502, 15, 26, "Input"],
Cell[51871, 1401, 312, 7, 42, "Text"],
Cell[52186, 1410, 194, 5, 24, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[52417, 1420, 26, 0, 54, "Section"],
Cell[52446, 1422, 905, 15, 135, "Text"]
}, Closed]]
}, Open  ]]
}
]
*)

(* End of internal cache information *)