<!--{{{-->
<link rel='alternate' type='application/rss+xml' title='RSS' href='index.xml' />
<!--}}}-->
Background: #fff
Foreground: #000
PrimaryPale: #8cf
PrimaryLight: #18f
PrimaryMid: #04b
PrimaryDark: #014
SecondaryPale: #ffc
SecondaryLight: #fe8
SecondaryMid: #db4
SecondaryDark: #841
TertiaryPale: #eee
TertiaryLight: #ccc
TertiaryMid: #999
TertiaryDark: #666
Error: #f88
/*{{{*/
body {background:[[ColorPalette::Background]]; color:[[ColorPalette::Foreground]];}

a {color:[[ColorPalette::PrimaryMid]];}
a:hover {background-color:[[ColorPalette::PrimaryMid]]; color:[[ColorPalette::Background]];}
a img {border:0;}

h1,h2,h3,h4,h5,h6 {color:[[ColorPalette::SecondaryDark]]; background:transparent;}
h1 {border-bottom:2px solid [[ColorPalette::TertiaryLight]];}
h2,h3 {border-bottom:1px solid [[ColorPalette::TertiaryLight]];}

.button {color:[[ColorPalette::PrimaryDark]]; border:1px solid [[ColorPalette::Background]];}
.button:hover {color:[[ColorPalette::PrimaryDark]]; background:[[ColorPalette::SecondaryLight]]; border-color:[[ColorPalette::SecondaryMid]];}
.button:active {color:[[ColorPalette::Background]]; background:[[ColorPalette::SecondaryMid]]; border:1px solid [[ColorPalette::SecondaryDark]];}

.header {background:[[ColorPalette::PrimaryMid]];}
.headerShadow {color:[[ColorPalette::Foreground]];}
.headerShadow a {font-weight:normal; color:[[ColorPalette::Foreground]];}
.headerForeground {color:[[ColorPalette::Background]];}
.headerForeground a {font-weight:normal; color:[[ColorPalette::PrimaryPale]];}

.tabSelected {color:[[ColorPalette::PrimaryDark]];
	background:[[ColorPalette::TertiaryPale]];
	border-left:1px solid [[ColorPalette::TertiaryLight]];
	border-top:1px solid [[ColorPalette::TertiaryLight]];
	border-right:1px solid [[ColorPalette::TertiaryLight]];
}
.tabUnselected {color:[[ColorPalette::Background]]; background:[[ColorPalette::TertiaryMid]];}
.tabContents {color:[[ColorPalette::PrimaryDark]]; background:[[ColorPalette::TertiaryPale]]; border:1px solid [[ColorPalette::TertiaryLight]];}
.tabContents .button {border:0;}

#sidebar {}
#sidebarOptions input {border:1px solid [[ColorPalette::PrimaryMid]];}
#sidebarOptions .sliderPanel {background:[[ColorPalette::PrimaryPale]];}
#sidebarOptions .sliderPanel a {border:none;color:[[ColorPalette::PrimaryMid]];}
#sidebarOptions .sliderPanel a:hover {color:[[ColorPalette::Background]]; background:[[ColorPalette::PrimaryMid]];}
#sidebarOptions .sliderPanel a:active {color:[[ColorPalette::PrimaryMid]]; background:[[ColorPalette::Background]];}

.wizard {background:[[ColorPalette::PrimaryPale]]; border:1px solid [[ColorPalette::PrimaryMid]];}
.wizard h1 {color:[[ColorPalette::PrimaryDark]]; border:none;}
.wizard h2 {color:[[ColorPalette::Foreground]]; border:none;}
.wizardStep {background:[[ColorPalette::Background]]; color:[[ColorPalette::Foreground]];
	border:1px solid [[ColorPalette::PrimaryMid]];}
.wizardStep.wizardStepDone {background:[[ColorPalette::TertiaryLight]];}
.wizardFooter {background:[[ColorPalette::PrimaryPale]];}
.wizardFooter .status {background:[[ColorPalette::PrimaryDark]]; color:[[ColorPalette::Background]];}
.wizard .button {color:[[ColorPalette::Foreground]]; background:[[ColorPalette::SecondaryLight]]; border: 1px solid;
	border-color:[[ColorPalette::SecondaryPale]] [[ColorPalette::SecondaryDark]] [[ColorPalette::SecondaryDark]] [[ColorPalette::SecondaryPale]];}
.wizard .button:hover {color:[[ColorPalette::Foreground]]; background:[[ColorPalette::Background]];}
.wizard .button:active {color:[[ColorPalette::Background]]; background:[[ColorPalette::Foreground]]; border: 1px solid;
	border-color:[[ColorPalette::PrimaryDark]] [[ColorPalette::PrimaryPale]] [[ColorPalette::PrimaryPale]] [[ColorPalette::PrimaryDark]];}

.wizard .notChanged {background:transparent;}
.wizard .changedLocally {background:#80ff80;}
.wizard .changedServer {background:#8080ff;}
.wizard .changedBoth {background:#ff8080;}
.wizard .notFound {background:#ffff80;}
.wizard .putToServer {background:#ff80ff;}
.wizard .gotFromServer {background:#80ffff;}

#messageArea {border:1px solid [[ColorPalette::SecondaryMid]]; background:[[ColorPalette::SecondaryLight]]; color:[[ColorPalette::Foreground]];}
#messageArea .button {color:[[ColorPalette::PrimaryMid]]; background:[[ColorPalette::SecondaryPale]]; border:none;}

.popupTiddler {background:[[ColorPalette::TertiaryPale]]; border:2px solid [[ColorPalette::TertiaryMid]];}

.popup {background:[[ColorPalette::TertiaryPale]]; color:[[ColorPalette::TertiaryDark]]; border-left:1px solid [[ColorPalette::TertiaryMid]]; border-top:1px solid [[ColorPalette::TertiaryMid]]; border-right:2px solid [[ColorPalette::TertiaryDark]]; border-bottom:2px solid [[ColorPalette::TertiaryDark]];}
.popup hr {color:[[ColorPalette::PrimaryDark]]; background:[[ColorPalette::PrimaryDark]]; border-bottom:1px;}
.popup li.disabled {color:[[ColorPalette::TertiaryMid]];}
.popup li a, .popup li a:visited {color:[[ColorPalette::Foreground]]; border: none;}
.popup li a:hover {background:[[ColorPalette::SecondaryLight]]; color:[[ColorPalette::Foreground]]; border: none;}
.popup li a:active {background:[[ColorPalette::SecondaryPale]]; color:[[ColorPalette::Foreground]]; border: none;}
.popupHighlight {background:[[ColorPalette::Background]]; color:[[ColorPalette::Foreground]];}
.listBreak div {border-bottom:1px solid [[ColorPalette::TertiaryDark]];}

.tiddler .defaultCommand {font-weight:bold;}

.shadow .title {color:[[ColorPalette::TertiaryDark]];}

.title {color:[[ColorPalette::SecondaryDark]];}
.subtitle {color:[[ColorPalette::TertiaryDark]];}

.toolbar {color:[[ColorPalette::PrimaryMid]];}
.toolbar a {color:[[ColorPalette::TertiaryLight]];}
.selected .toolbar a {color:[[ColorPalette::TertiaryMid]];}
.selected .toolbar a:hover {color:[[ColorPalette::Foreground]];}

.tagging, .tagged {border:1px solid [[ColorPalette::TertiaryPale]]; background-color:[[ColorPalette::TertiaryPale]];}
.selected .tagging, .selected .tagged {background-color:[[ColorPalette::TertiaryLight]]; border:1px solid [[ColorPalette::TertiaryMid]];}
.tagging .listTitle, .tagged .listTitle {color:[[ColorPalette::PrimaryDark]];}
.tagging .button, .tagged .button {border:none;}

.footer {color:[[ColorPalette::TertiaryLight]];}
.selected .footer {color:[[ColorPalette::TertiaryMid]];}

.error, .errorButton {color:[[ColorPalette::Foreground]]; background:[[ColorPalette::Error]];}
.warning {color:[[ColorPalette::Foreground]]; background:[[ColorPalette::SecondaryPale]];}
.lowlight {background:[[ColorPalette::TertiaryLight]];}

.zoomer {background:none; color:[[ColorPalette::TertiaryMid]]; border:3px solid [[ColorPalette::TertiaryMid]];}

.imageLink, #displayArea .imageLink {background:transparent;}

.annotation {background:[[ColorPalette::SecondaryLight]]; color:[[ColorPalette::Foreground]]; border:2px solid [[ColorPalette::SecondaryMid]];}

.viewer .listTitle {list-style-type:none; margin-left:-2em;}
.viewer .button {border:1px solid [[ColorPalette::SecondaryMid]];}
.viewer blockquote {border-left:3px solid [[ColorPalette::TertiaryDark]];}

.viewer table, table.twtable {border:2px solid [[ColorPalette::TertiaryDark]];}
.viewer th, .viewer thead td, .twtable th, .twtable thead td {background:[[ColorPalette::SecondaryMid]]; border:1px solid [[ColorPalette::TertiaryDark]]; color:[[ColorPalette::Background]];}
.viewer td, .viewer tr, .twtable td, .twtable tr {border:1px solid [[ColorPalette::TertiaryDark]];}

.viewer pre {border:1px solid [[ColorPalette::SecondaryLight]]; background:[[ColorPalette::SecondaryPale]];}
.viewer code {color:[[ColorPalette::SecondaryDark]];}
.viewer hr {border:0; border-top:dashed 1px [[ColorPalette::TertiaryDark]]; color:[[ColorPalette::TertiaryDark]];}

.highlight, .marked {background:[[ColorPalette::SecondaryLight]];}

.editor input {border:1px solid [[ColorPalette::PrimaryMid]];}
.editor textarea {border:1px solid [[ColorPalette::PrimaryMid]]; width:100%;}
.editorFooter {color:[[ColorPalette::TertiaryMid]];}
.readOnly {background:[[ColorPalette::TertiaryPale]];}

#backstageArea {background:[[ColorPalette::Foreground]]; color:[[ColorPalette::TertiaryMid]];}
#backstageArea a {background:[[ColorPalette::Foreground]]; color:[[ColorPalette::Background]]; border:none;}
#backstageArea a:hover {background:[[ColorPalette::SecondaryLight]]; color:[[ColorPalette::Foreground]]; }
#backstageArea a.backstageSelTab {background:[[ColorPalette::Background]]; color:[[ColorPalette::Foreground]];}
#backstageButton a {background:none; color:[[ColorPalette::Background]]; border:none;}
#backstageButton a:hover {background:[[ColorPalette::Foreground]]; color:[[ColorPalette::Background]]; border:none;}
#backstagePanel {background:[[ColorPalette::Background]]; border-color: [[ColorPalette::Background]] [[ColorPalette::TertiaryDark]] [[ColorPalette::TertiaryDark]] [[ColorPalette::TertiaryDark]];}
.backstagePanelFooter .button {border:none; color:[[ColorPalette::Background]];}
.backstagePanelFooter .button:hover {color:[[ColorPalette::Foreground]];}
#backstageCloak {background:[[ColorPalette::Foreground]]; opacity:0.6; filter:alpha(opacity=60);}
/*}}}*/
/*{{{*/
* html .tiddler {height:1%;}

body {font-size:.75em; font-family:arial,helvetica; margin:0; padding:0;}

h1,h2,h3,h4,h5,h6 {font-weight:bold; text-decoration:none;}
h1,h2,h3 {padding-bottom:1px; margin-top:1.2em;margin-bottom:0.3em;}
h4,h5,h6 {margin-top:1em;}
h1 {font-size:1.35em;}
h2 {font-size:1.25em;}
h3 {font-size:1.1em;}
h4 {font-size:1em;}
h5 {font-size:.9em;}

hr {height:1px;}

a {text-decoration:none;}

dt {font-weight:bold;}

ol {list-style-type:decimal;}
ol ol {list-style-type:lower-alpha;}
ol ol ol {list-style-type:lower-roman;}
ol ol ol ol {list-style-type:decimal;}
ol ol ol ol ol {list-style-type:lower-alpha;}
ol ol ol ol ol ol {list-style-type:lower-roman;}
ol ol ol ol ol ol ol {list-style-type:decimal;}

.txtOptionInput {width:11em;}

#contentWrapper .chkOptionInput {border:0;}

.externalLink {text-decoration:underline;}

.indent {margin-left:3em;}
.outdent {margin-left:3em; text-indent:-3em;}
code.escaped {white-space:nowrap;}

.tiddlyLinkExisting {font-weight:bold;}
.tiddlyLinkNonExisting {font-style:italic;}

/* the 'a' is required for IE, otherwise it renders the whole tiddler in bold */
a.tiddlyLinkNonExisting.shadow {font-weight:bold;}

#mainMenu .tiddlyLinkExisting,
	#mainMenu .tiddlyLinkNonExisting,
	#sidebarTabs .tiddlyLinkNonExisting {font-weight:normal; font-style:normal;}
#sidebarTabs .tiddlyLinkExisting {font-weight:bold; font-style:normal;}

.header {position:relative;}
.header a:hover {background:transparent;}
.headerShadow {position:relative; padding:4.5em 0 1em 1em; left:-1px; top:-1px;}
.headerForeground {position:absolute; padding:4.5em 0 1em 1em; left:0; top:0;}

.siteTitle {font-size:3em;}
.siteSubtitle {font-size:1.2em;}

#mainMenu {position:absolute; left:0; width:10em; text-align:right; line-height:1.6em; padding:1.5em 0.5em 0.5em 0.5em; font-size:1.1em;}

#sidebar {position:absolute; right:3px; width:16em; font-size:.9em;}
#sidebarOptions {padding-top:0.3em;}
#sidebarOptions a {margin:0 0.2em; padding:0.2em 0.3em; display:block;}
#sidebarOptions input {margin:0.4em 0.5em;}
#sidebarOptions .sliderPanel {margin-left:1em; padding:0.5em; font-size:.85em;}
#sidebarOptions .sliderPanel a {font-weight:bold; display:inline; padding:0;}
#sidebarOptions .sliderPanel input {margin:0 0 0.3em 0;}
#sidebarTabs .tabContents {width:15em; overflow:hidden;}

.wizard {padding:0.1em 1em 0 2em;}
.wizard h1 {font-size:2em; font-weight:bold; background:none; padding:0; margin:0.4em 0 0.2em;}
.wizard h2 {font-size:1.2em; font-weight:bold; background:none; padding:0; margin:0.4em 0 0.2em;}
.wizardStep {padding:1em 1em 1em 1em;}
.wizard .button {margin:0.5em 0 0; font-size:1.2em;}
.wizardFooter {padding:0.8em 0.4em 0.8em 0;}
.wizardFooter .status {padding:0 0.4em; margin-left:1em;}
.wizard .button {padding:0.1em 0.2em;}

#messageArea {position:fixed; top:2em; right:0; margin:0.5em; padding:0.5em; z-index:2000; _position:absolute;}
.messageToolbar {display:block; text-align:right; padding:0.2em;}
#messageArea a {text-decoration:underline;}

.tiddlerPopupButton {padding:0.2em;}
.popupTiddler {position: absolute; z-index:300; padding:1em; margin:0;}

.popup {position:absolute; z-index:300; font-size:.9em; padding:0; list-style:none; margin:0;}
.popup .popupMessage {padding:0.4em;}
.popup hr {display:block; height:1px; width:auto; padding:0; margin:0.2em 0;}
.popup li.disabled {padding:0.4em;}
.popup li a {display:block; padding:0.4em; font-weight:normal; cursor:pointer;}
.listBreak {font-size:1px; line-height:1px;}
.listBreak div {margin:2px 0;}

.tabset {padding:1em 0 0 0.5em;}
.tab {margin:0 0 0 0.25em; padding:2px;}
.tabContents {padding:0.5em;}
.tabContents ul, .tabContents ol {margin:0; padding:0;}
.txtMainTab .tabContents li {list-style:none;}
.tabContents li.listLink { margin-left:.75em;}

#contentWrapper {display:block;}
#splashScreen {display:none;}

#displayArea {margin:1em 17em 0 14em;}

.toolbar {text-align:right; font-size:.9em;}

.tiddler {padding:1em 1em 0;}

.missing .viewer,.missing .title {font-style:italic;}

.title {font-size:1.6em; font-weight:bold;}

.missing .subtitle {display:none;}
.subtitle {font-size:1.1em;}

.tiddler .button {padding:0.2em 0.4em;}

.tagging {margin:0.5em 0.5em 0.5em 0; float:left; display:none;}
.isTag .tagging {display:block;}
.tagged {margin:0.5em; float:right;}
.tagging, .tagged {font-size:0.9em; padding:0.25em;}
.tagging ul, .tagged ul {list-style:none; margin:0.25em; padding:0;}
.tagClear {clear:both;}

.footer {font-size:.9em;}
.footer li {display:inline;}

.annotation {padding:0.5em; margin:0.5em;}

* html .viewer pre {width:99%; padding:0 0 1em 0;}
.viewer {line-height:1.4em; padding-top:0.5em;}
.viewer .button {margin:0 0.25em; padding:0 0.25em;}
.viewer blockquote {line-height:1.5em; padding-left:0.8em;margin-left:2.5em;}
.viewer ul, .viewer ol {margin-left:0.5em; padding-left:1.5em;}

.viewer table, table.twtable {border-collapse:collapse; margin:0.8em 1.0em;}
.viewer th, .viewer td, .viewer tr,.viewer caption,.twtable th, .twtable td, .twtable tr,.twtable caption {padding:3px;}
table.listView {font-size:0.85em; margin:0.8em 1.0em;}
table.listView th, table.listView td, table.listView tr {padding:0 3px 0 3px;}

.viewer pre {padding:0.5em; margin-left:0.5em; font-size:1.2em; line-height:1.4em; overflow:auto;}
.viewer code {font-size:1.2em; line-height:1.4em;}

.editor {font-size:1.1em;}
.editor input, .editor textarea {display:block; width:100%; font:inherit;}
.editorFooter {padding:0.25em 0; font-size:.9em;}
.editorFooter .button {padding-top:0; padding-bottom:0;}

.fieldsetFix {border:0; padding:0; margin:1px 0px;}

.zoomer {font-size:1.1em; position:absolute; overflow:hidden;}
.zoomer div {padding:1em;}

* html #backstage {width:99%;}
* html #backstageArea {width:99%;}
#backstageArea {display:none; position:relative; overflow: hidden; z-index:150; padding:0.3em 0.5em;}
#backstageToolbar {position:relative;}
#backstageArea a {font-weight:bold; margin-left:0.5em; padding:0.3em 0.5em;}
#backstageButton {display:none; position:absolute; z-index:175; top:0; right:0;}
#backstageButton a {padding:0.1em 0.4em; margin:0.1em;}
#backstage {position:relative; width:100%; z-index:50;}
#backstagePanel {display:none; z-index:100; position:absolute; width:90%; margin-left:3em; padding:1em;}
.backstagePanelFooter {padding-top:0.2em; float:right;}
.backstagePanelFooter a {padding:0.2em 0.4em;}
#backstageCloak {display:none; z-index:20; position:absolute; width:100%; height:100px;}

.whenBackstage {display:none;}
.backstageVisible .whenBackstage {display:block;}
/*}}}*/
/***
StyleSheet for use when a translation requires any css style changes.
This StyleSheet can be used directly by languages such as Chinese, Japanese and Korean which need larger font sizes.
***/
/*{{{*/
body {font-size:0.8em;}
#sidebarOptions {font-size:1.05em;}
#sidebarOptions a {font-style:normal;}
#sidebarOptions .sliderPanel {font-size:0.95em;}
.subtitle {font-size:0.8em;}
.viewer table.listView {font-size:0.95em;}
/*}}}*/
/*{{{*/
@media print {
#mainMenu, #sidebar, #messageArea, .toolbar, #backstageButton, #backstageArea {display: none !important;}
#displayArea {margin: 1em 1em 0em;}
noscript {display:none;} /* Fixes a feature in Firefox 1.5.0.2 where print preview displays the noscript content */
}
/*}}}*/
<!--{{{-->
<div class='header' macro='gradient vert [[ColorPalette::PrimaryLight]] [[ColorPalette::PrimaryMid]]'>
<div class='headerShadow'>
<span class='siteTitle' refresh='content' tiddler='SiteTitle'></span>&nbsp;
<span class='siteSubtitle' refresh='content' tiddler='SiteSubtitle'></span>
</div>
<div class='headerForeground'>
<span class='siteTitle' refresh='content' tiddler='SiteTitle'></span>&nbsp;
<span class='siteSubtitle' refresh='content' tiddler='SiteSubtitle'></span>
</div>
</div>
<div id='mainMenu' refresh='content' tiddler='MainMenu'></div>
<div id='sidebar'>
<div id='sidebarOptions' refresh='content' tiddler='SideBarOptions'></div>
<div id='sidebarTabs' refresh='content' force='true' tiddler='SideBarTabs'></div>
</div>
<div id='displayArea'>
<div id='messageArea'></div>
<div id='tiddlerDisplay'></div>
</div>
<!--}}}-->
<!--{{{-->
<div class='toolbar' macro='toolbar [[ToolbarCommands::ViewToolbar]]'></div>
<div class='title' macro='view title'></div>
<div class='subtitle'><span macro='view modifier link'></span>, <span macro='view modified date'></span> (<span macro='message views.wikified.createdPrompt'></span> <span macro='view created date'></span>)</div>
<div class='tagging' macro='tagging'></div>
<div class='tagged' macro='tags'></div>
<div class='viewer' macro='view text wikified'></div>
<div class='tagClear'></div>
<!--}}}-->
<!--{{{-->
<div class='toolbar' macro='toolbar [[ToolbarCommands::EditToolbar]]'></div>
<div class='title' macro='view title'></div>
<div class='editor' macro='edit title'></div>
<div macro='annotations'></div>
<div class='editor' macro='edit text'></div>
<div class='editor' macro='edit tags'></div><div class='editorFooter'><span macro='message views.editor.tagPrompt'></span><span macro='tagChooser excludeLists'></span></div>
<!--}}}-->
To get started with this blank [[TiddlyWiki]], you'll need to modify the following tiddlers:
* [[SiteTitle]] & [[SiteSubtitle]]: The title and subtitle of the site, as shown above (after saving, they will also appear in the browser title bar)
* [[MainMenu]]: The menu (usually on the left)
* [[DefaultTiddlers]]: Contains the names of the tiddlers that you want to appear when the TiddlyWiki is opened
You'll also need to enter your username for signing your edits: <<option txtUserName>>
These [[InterfaceOptions]] for customising [[TiddlyWiki]] are saved in your browser

Your username for signing your edits. Write it as a [[WikiWord]] (eg [[JoeBloggs]])

<<option txtUserName>>
<<option chkSaveBackups>> [[SaveBackups]]
<<option chkAutoSave>> [[AutoSave]]
<<option chkRegExpSearch>> [[RegExpSearch]]
<<option chkCaseSensitiveSearch>> [[CaseSensitiveSearch]]
<<option chkAnimate>> [[EnableAnimations]]

----
Also see [[AdvancedOptions]]
<<importTiddlers>>
!! Literature
# Wm. G. Hoover and C.G. Hoover. Ergodicity of a ~Time-Reversibly Thermostated Harmonic Oscillator and the 2014 Ian Snook Prize. [[(2014)|http://arxiv.org/abs/1408.0256]]. Shuichi Nose opened up a new world of atomistic simulation in 1984. He formulated a Hamiltonian tailored to generate Gibbs' canonical distribution dynamically. This clever idea bridged the gap between microcanonical molecular dynamics and canonical statistical mechanics.
* The status quo for learning C++: http://www.learncpp.com/ 
* Great option for those who are confident with other programming languages but need to transfer that to knowledge C++: http://www.cplusplus.com/doc/tutorial/
* https://www.sololearn.com/Course/CPlusPlus/  uses a very interactive learning structure with short readings and examples follows by questions which one must get right to move to the next part of the tutorial.
* http://www.cprogramming.com/tutorial/c++-tutorial.html. 
* http://www.cprogramming.com/tutorial/c-tutorial.html
* http://www.learn-c.org/
# [[Basic introduction|http://homepages.spa.umn.edu/~vinals/classes/phys4041/phys4041_notes.pdf]] to scientific programming elements which we will be using in class.
# [[Computer Hardware]].
# [[Programming Tips]].
# Text of [[Numerical Recipes|http://www.aip.de/groups/soe/local/numres/]].
# C [[code|https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_sin.c;hb=HEAD#l281]] for computing sin(x) or cos(x).
# [[Central Limit Theorem|http://www.math.csusb.edu/faculty/stanton/probstat/clt.html]]. An applet demonstrating the limit.
# LU Decomposition (John Burkardt, Florida State University, http://people.sc.fsu.edu/~jburkardt/classes/gateway_2014/gateway_2014.html ). [[pdf|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/LUDecomposition.pdf]].
# The Conjugate Gradient [[method|http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf]].
# Examples of PDE solutions (Java applets).
## [[Poisson equation|http://www.um.es/fem/EjsWiki/Main/ExamplesPoisson]].
## [[One dimensional, time dependent Schrodinger equation|http://www.math.stonybrook.edu/~ccc/WaveSim.html]].
## [[Assorted examples|http://www.falstad.com/mathphysics.html]]. 2D waves and Poisson.
## [[One dimensional diffusion equation|https://ocw.mit.edu/ans7870/18/18.03/s06/tools/HeatEquation.html]].
## [[A collection of links|http://www.math.lamar.edu/faculty/maesumi/PDE1.html]].
# [[Ewald Sums|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/EwaldSum.pdf]].
# Molecular Dynamics simulations.
## [[An introduction|http://www.youtube.com/watch?v=lLFEqKl3sm4&list=PLY5jumQrwEXXlVEg1bsUR72yc8v3t4dJc]]
## [[Water|http://www.youtube.com/watch?v=x8Atqz5YvzQ]].
## [[Peptide|http://www.youtube.com/watch?v=WScPbvUwDno]]
## [[Ice melting|https://www.youtube.com/watch?v=NQhjAtCKghE]].
## A simple deterministic and time reversal invariant thermostat. Henk van [[Beijeren|http://arxiv.org/abs/1411.2983]].
# [[Principal Component Analysis example]].
# [[K-means]].
# [[Expectation Maximization Clustering]] example.
# [[Example of pseudo-code|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2016/p_code.pdf]].
The innards of a computer begin in its motherboard. This is the base where the CPU, memory, chipset, network, and connectivity slots are located. An annotated image can be found [[here|http://i.nextmedia.com.au/features/i7-board-diagram.jpg]].

As shown in the motherboard images, the key part of a computer is the so called processor. They are built in large microfabriaction facilities, contain approximately 1 billion transistors, with individual features down to about 20 nm (and shrinking, 7 nm is the latest technology).

A relativevely new modern architecture is the Intel Sandy Bridge processor (with 32 nm features). The processor containes 6 cores or processing elements that share the same 'die", its own cache and memory controleler. An image can be found at http://www.hardwarezone.com.sg/product-intel-core-i7-3960x-extreme-edition

More details can be found in: [[A beginner's guide to High Performance Computing|http://www.shodor.org/media/content/petascale/materials/UPModules/beginnersGuideHPC/moduleDocument_pdf.pdf]] by Rubin Landau, Oregon State U.
[[Home]]
Given the data $\vec{x} = \{1.0,1.3,2.2,2.6,2.8,5.0,7.3,7.4,7.5,7.7,7.9\}$, two Gaussian clusters $k=2$, initial guesses for cluster means $\mu_{1} = 6.63$ and $\mu_{2} = 7.57$, guesses for variances $\sigma_{1} = \sigma_{2} = 1.0$ and weights $c_{1} = c_{2} = 0.5$, the following code performs the iteration:
{{{
PROGRAM em_clustering

! Example 13.4, p. 346. Data mining and analysis. Mohamed Zaki and Wagner 
! Meira. Cambridge University Press, 2014.

REAL(8) :: x(11),mu(2),sigma(2),c(2)
REAL(8) :: wg(2,11),w(2,11)

DATA x / 1.0,1.3,2.2,2.6,2.8,5.0,7.3,7.4,7.5,7.7,7.9 /

tpi2=sqrt(8.0*atan(1.0))

mu(1)=6.63
mu(2)=7.57
sigma(1)=1.0
sigma(2)=1.0
c(1)=0.5
c(2)=0.5

! Set up the cluster Gaussians for the current values of the parameters
DO i=1,2
DO j=1,11
wg(i,j)=exp(-(x(j)-mu(i))**2/(2.0*sigma(i)**2))/(tpi2*sigma(i))
END DO
END DO

! Use Bayes theorem to compute the probabilities of the parameters given the data
DO i=1,2
DO j=1,11
w(i,j)=wg(i,j)*c(i)/(wg(1,j)*c(1)+wg(2,j)*c(2))
END DO
END DO

! Recompute means, variances and weights with the newly computed probailities of the parameters
DO i=1,2
sumw=0.0
mu(i)=0.0
DO j=1,11
sumw=sumw+w(i,j)
mu(i)=mu(i)+w(i,j)*x(j)
END DO
mu(i)=mu(i)/sumw
END DO

DO i=1,2
sumw=0.0
sigma(i)=0.0
DO j=1,11
sumw=sumw+w(i,j)
sigma(i)=sigma(i)+w(i,j)*(x(j)-mu(i))**2
END DO
sigma(i)=sigma(i)/sumw
END DO

DO i=1,2
c(i)=0.0
DO j=1,11
c(i)=c(i)+w(i,j)
END DO
c(i)=c(i)/11.0
END DO

DO i=1,2
write(6,9999) i, mu(i),sigma(i),c(i)
9999 FORMAT(i4,4x,3f10.3)
END DO

END
}}}

[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/em_clustering.jpg]]

[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/em_clustering_2d.jpg]]
[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/images/PANO_20130902_180940_rotated2.jpg]]

!General Information
* [[Computing resources in the School of Physics|http://zzz.physics.umn.edu/computing/home]].
* Data plotting and visualization: [[gnuplot|http://www.gnuplot.info/]] and [[xmgrace|http://plasma-gate.weizmann.ac.il/Grace/doc/UsersGuide.html]]. Both are available on School computers and the CSE Labs.
* MATLAB can also be used for data manipulation and visualization (although do not use MATLAB for your homework assignments). The University has a site license for the package ([[information|http://zzz.physics.umn.edu/computing/software/matlab]]). To use, just type matlab in any School computer.
* [[Extended course description|http://www.aem.umn.edu/~shield/itcc/minutes/2011-4-19/phys4041.html]].
! Programming language introductions
To follow current standards in scientific computing, you must use at least one of the following high level languages in homework assignments or lab sessions: C, C++ or FORTRAN. You are welcome to use MATLAB, Python, Java, etc. for your own development needs, but only the three given languages are acceptable in the class.
* [[FORTRAN tutorials|http://fortranwiki.org/fortran/show/Tutorials]]. FORTRAN reference [[here|http://www.am.qub.ac.uk/users/d.dundas/msci_workshop/fortran/notes/fortran90.pdf]] and [[here|http://www.personal.psu.edu/jhm/f90/lectures/]]. See also [[AMath 483|http://faculty.washington.edu/rjl/uwamath583s11/sphinx/notes/html/]] at University of Washington, including an introduction to Unix.
* [[C tutorials|CTutorials]], [[C++ tutorials|C++Tutorials]].
* Analysis of [[common features|http://www.featflow.de/beta/en/software/featflow2/tutorial/tutorial_lang.html]] between FORTRAN and C++.
* [[Python|http://faculty.washington.edu/rjl/uwamath583s11/sphinx/notes/html/]].
* [[Introduction|http://www.cs.rpi.edu/~szymansk/OOF90/F90_Objects.html]] to object oriented FORTRAN 90.
* Intel Math Kernel Lirbrary, Reference Manual: https://software.intel.com/en-us/mkl_11.2_ref 
! CSE Computer Lab information
* General [[information|http://cselabs.umn.edu/]], including [[account creation|http://cselabs.umn.edu/accounts]].
* Local access to Unix environment from Keller 1-200.
** At the Start menu (left hand side bottom), start the ~FastX 2 application.
** Click on + sign (right side of menu), and click on web.
** Name: enter any name for session. Host: https://csel-vole-nn.cselabs.umn.edu:3443 nn are the two digits shown on the label of your screen monitor (lower left). User: ~YourUsername. Save.
** Double click on the connection just created. Enter your username and password at the prompts. Click continue on the next window that appears. 
** Click on + sign on top right (create a connection). Eventually a "deskpace image" will appear. Double click on it. Now you have a Unix session running. You can maximize the window if you want it to take the whole screen.
** To terminate, just click on window icon x
** At the bottom of the Unix desktop, you will find an icon for a "terminal" and one for a browser. You can click on them to open.
* Remote access to Labs: 
** Browser based access to Unix environment: http://vole.cselabs.umn.edu
** General from Unix shell: ssh -X vole.cse.umn.edu -l your_~CSELabs_username ([[further information on ssh|http://help.cselabs.umn.edu/offsite/ssh]]). Also [[list| http://cselabs.umn.edu/labs/unix_machines]] of hosts available.
** [[Further information|http://help.cselabs.umn.edu/offsite]] about methods for offsite access.
* ''FORTRAN'' (type man gfortran for information and options)
** gfortran filename.f90
* ''C'' (type man g++ or man gcc for information and options)
** g++ filename.c or gcc filename.c -lm
! ~GitHub
The version control and repository environment ~GitHub will be used to make code and other class materials available to the class.
* Access: https://github.umn.edu .
* Log in and email the instructor when you have done so. The instructor will then authorize you to access the repository for this class.
* Initialize. Go to (or create) the subdirectory in the CSE Labs computer that you want to use for this class and type 
{{{
git clone https://github.umn.edu/jvinals/PHYS4041-2018.git
}}}
You can also clone (make a copy) of all the repository content in your home computer or laptop by issuing the same command. Keep in mind that you will need to have git installed on that computer. Versions of git exist for all modern operating systems, including phones (if you care about such things).

To update the local contents on your directory on the CSE Labs host, go to that directory and type
{{{
git pull
}}}
{{justifyright{
''Due Tuesday, September 16, 2014''
}}}

@@color(blue):Important note:@@ Code needs to be submitted by email to the TA. Other material: Hardcopies with plots, calculations, etc. (but ''not code'') are to be submitted directly to the instructor. It is also possible to submit graphs in image form or scripts to generate them directly to the TA in electronic form. It is best if you group all files into a single archive and email the archive. For example, create a directory homework1, copy all files that need to be submitted to it and create the archive by typing:
* % mkdir homework1
* % cp [filename1] [filename2] ... [filenamen] homework1/
* % tar cvf homework1.tar homework1
and email file homework1.tar

<hr>

# Write a subroutine or a function that takes an array of real numbers as argument, and that returns +1 if the sum of the array elements is positive, 0 if the sum of array elements is zero, and -1 if the sum is negative. Write a main program that uses this subroutine or function on a particular example of your choice.
# Consider the following system of two equations with two unknowns \begin{eqnarray}
3x + ay & = & 10, \\
5x + by & = & 20,
\end{eqnarray}
with $a = 2.100 \pm 5 \times 10^{-4}$ and $b = 3.300 \pm 5 \times 10^{-4}$. Solve for $y$ and find the error in $y$.
# Given a sequence of real numbers $\{x_{1}, x_{2}, ... x_{n} \}$, devise and describe an iterative algorithm to sort the sequence in (a) ascending value, (b) ascending absolute value. Write a computer program to implement your algorithm, and run it on a test sequence of your choosing.
# Write some code to establish the number of digits of precision available in both single and double precision floating point arithmetic in the computer that you are using. To accomplish this task, obtain the following results with successively smaller values of $\epsilon$: $x = 1 + \epsilon$, $y = 1 - \epsilon$ and $z = x - y$. Mathematically, $z = 2\epsilon$, but as $\epsilon$ becomes sufficiently small, the code will return $z = 0$ instead. Find for both cases of single and double precision how small can $\epsilon$ be until you obtain $z = 0$; this is the precision of the computer you are using.
# Write some code that uses single precision arithmetic to evaluate the sum 
$$ S_{N} = \sum_{n=1}^{2N} (-1)^{n} \frac{n}{n+1} = - \sum_{n=1}^{N} \frac{2n-1}{2n} + \sum_{n=1}^{N} \frac{2n}{2n+1} = \sum_{n=1}^{N} \frac{1}{2n(2n+1)}.
$$
## Evaluate the sum $S_{N}$ by using the three algebraic expressions given above when N = 100, 1000, 10000. Comment on any differences you may find among the three calculations.
## Using the third algebraic formula determine how many terms would you require in the calculation to make sure that the sum has converged to an accuracy of $10^{-7}$.

[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/PHYS4041_Homework1.pdf]]
{{justifyright{
''Due Friday, September 23, 2016''
}}}

@@color(blue):Important note:@@ Code needs to be submitted by email to the TA. Other material: Hardcopies with plots, calculations, etc. (but ''not code'') are to be submitted directly to the instructor or TA. It is also possible to submit graphs directly to the TA in electronic form. It is best if you group all files into a single archive and email the archive. For example, create a directory homework1, copy all files that need to be submitted to it and create the archive by typing:
{{{
mkdir homework1
cp [filename1] [filename2] ... [filenamen] homework1/
tar cvf homework1.tar homework1
}}}
and email file homework1.tar

<hr>

# Write a subroutine or a function that takes an array of real numbers as argument, and that returns +1 if the sum of the array elements is positive, 0 if the sum of array elements is zero, and -1 if the sum is negative. Write a main program that uses this subroutine or function on a particular example of your choice.
# Consider the following system of two equations with two unknowns \begin{eqnarray}
3x + ay & = & 10, \\
5x + by & = & 20,
\end{eqnarray}
with $a = 2.100 \pm 5 \times 10^{-4}$ and $b = 3.300 \pm 5 \times 10^{-4}$. Solve for $y$ and find the error in $y$.
# Given a sequence of real numbers $\{x_{1}, x_{2}, ... x_{n} \}$, devise and describe an iterative algorithm to sort the sequence in (a) ascending value, (b) ascending absolute value. Write a computer program to implement your algorithm, and run it on a test sequence of your choosing.
# Write some code to establish the number of digits of precision available in both single and double precision floating point arithmetic in the computer that you are using. To accomplish this task, obtain the following results with successively smaller values of $\epsilon$: $x = 1 + \epsilon$, $y = 1 - \epsilon$ and $z = x - y$. Mathematically, $z = 2\epsilon$, but as $\epsilon$ becomes sufficiently small, the code will return $z = 0$ instead. Find for both cases of single and double precision how small can $\epsilon$ be until you obtain $z = 0$; this is the precision of the computer you are using.
# Write some code that uses single precision arithmetic to evaluate the sum 
$$ S_{N} = \sum_{n=1}^{2N} (-1)^{n} \frac{n}{n+1} = - \sum_{n=1}^{N} \frac{2n-1}{2n} + \sum_{n=1}^{N} \frac{2n}{2n+1} = \sum_{n=1}^{N} \frac{1}{2n(2n+1)}.
$$
## Evaluate the sum $S_{N}$ by using the three algebraic expressions given above when N = 100, 1000, 10000. Comment on any differences you may find among the three calculations.
## Using the third algebraic formula determine how many terms would you require in the calculation to make sure that the sum has converged to an accuracy of $10^{-7}$.

[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2016/PHYS4041_Homework1.pdf]]
{{justifyright{
''Due Tuesday, September 18, 2018''
}}}

@@color(blue):Important note:@@ Code needs to be submitted by email to the TA. Other material: Hardcopies with plots, calculations, etc. (but ''not code'') are to be submitted directly to the instructor or TA. It is also possible to submit graphs directly to the TA in electronic form. It is best if you group all files into a single archive and email the archive. For example, create a directory homework1, copy all files that need to be submitted to it and create the archive by typing:
{{{
mkdir homework1
cp [filename1] [filename2] ... [filenamen] homework1/
tar cvf homework1.tar homework1
}}}
and email file homework1.tar

<hr>

!!! Question 1. Attainable accuracy in single and double precision arithmetic
Write and use code to establish the number of available digits in both single and double precision floating point arithmetic in the computer that you are using (state what computer it is). In order to do this, write a code that computes the following quantities as a function of $\epsilon$: $x= 1.0 + \epsilon$, $y = 1.0 - \epsilon$, $z = x -y$. The exact result is $z = 2 \epsilon$. However, as you run your code for smaller and smaller values of $\epsilon$, it will eventually return $z = 0$. At the point, the value of $\epsilon$ used is the precision of your computer.

!!! Question 2. Practice function coding
Write a subroutine or a function that takes an array of real numbers as argument, and that returns +1 if the sum of the array elements is positive, 0 if the sum of array elements is zero, and -1 if the sum is negative. Write a main program that uses this subroutine or function on a particular example of your choice.

!!! Question 3. Error propagation
Consider the following system of two equations with two unknowns \begin{eqnarray}
3x + ay & = & 10, \\
5x + by & = & 20,
\end{eqnarray}
with $a = 2.100 \pm 5 \times 10^{-4}$ and $b = 3.300 \pm 5 \times 10^{-4}$. Solve for $y$ and find the error in $y$.

!!! Question 4. Error accumulation in series evaluation
Write some code that uses single precision arithmetic to evaluate the sum 
$$ S_{N} = \sum_{n=1}^{2N} (-1)^{n} \frac{n}{n+1} = - \sum_{n=1}^{N} \frac{2n-1}{2n} + \sum_{n=1}^{N} \frac{2n}{2n+1} = \sum_{n=1}^{N} \frac{1}{2n(2n+1)}.
$$
## Evaluate the sum $S_{N}$ by using the three algebraic expressions given above when N = 100, 1000, 10000. Comment on any differences you may find among the three calculations.
## Using the third algebraic formula determine how many terms would you require in the calculation to make sure that the sum has converged to an accuracy of $10^{-7}$.

!!! Question 5. Algorithm to conduct an iteration
Given a sequence of real numbers $\{x_{1}, x_{2}, ... x_{n} \}$, devise and describe an iterative algorithm to sort the sequence in (a) ascending value, (b) ascending absolute value. Write a computer program to implement your algorithm, and run it on a test sequence of your choosing.
{{justifyright{
''Due Tuesday, September 23, 2014''
}}}

!!!!1. Bessel function approximation
Bessel functions are trancendental functions that can be obtained by summing a series. However, there are recursion formulae that make the computation more efficient, although, as seen in class, recursion can introduce its own numerical problems. Use the following recursion formulae to compute the Bessel function $J_{l}$ with $l = 0, 1, 2, 3, 4$ and plot the results for $x \in (0,12)$. The recursion formulae ("upwards" or "downwards") are \begin{eqnarray}
J_{l+1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l-1}(x) \\
J_{l-1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l+1}(x) \\
J_{0}(x) = \frac{\sin x}{x} & \; & J_{1}(x) = \frac{\sin x - x \cos x}{x^{2}}
\end{eqnarray}  Using the following table to help evaluate the accuracy of your approximations,
|! x |! $J_{3}$ |! $J_{5}$ |! $J_{8}$ |
| 0.10 | 9. 518519719 10$^{-6}$ | 9.616310231 10$^{-10}$ | 2.901200102 10$^{-16}$ |
| 1.00 | 9.006581118 10$^{-3}$ | 9.256115862 10$^{-5}$ | 2.826498802 10$^{-8}$ |
| 10.0 | -3.949584498 10$^{-2}$ | -5.553451162 10$^{-1}$ | 1.255780236 10$^{-1}$ |
** Compute the Bessel functions required by using the upward recursion.
** Compare your solution to the table. Argue how errors might have propagated and what effect this may have on the accuracy of the solution.
** Chose an arbitrary starting function $J_{L}(x)$ with large $L$ and use the iteration downward. Compare the results to the upward iteration for the Bessel functions asked above.

!!!!2. Spline fitting
We have found experimentally the pressure $p$ as a function of the molar volume ($v$) of a binary mixture at constant temperature and composition (in some arbitrary units). The results are shown in the following table.
|v | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 |
| p | 4.7 | 3.5 | 3.0 | 2.7 | 3.0 | 2.4 |
## Write your own program to find the cubic spline interpolation to $p=p(v)$ in this interval. The description of the algorithm can be found at http://mathworld.wolfram.com/CubicSpline.html
## Use the spline interpolation to find the value of the volume at which $p=3.25$.
## Accomplish the same task but download publicly available code on the web (scroll to the bottom of the page):
** Fortran version: http://people.sc.fsu.edu/~jburkardt/f_src/spline/spline.html
** C version: http://people.sc.fsu.edu/~jburkardt/c_src/spline/spline.html
** C++ version: http://people.sc.fsu.edu/~jburkardt/cpp_src/spline/spline.html

[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/PHYS4041_Homework2.pdf]]
{{justifyright{
''Due Friday, September 30, 2016''
}}}
!!! 1. Solution of simple linear system of equations (LAPACK)
Write a code (and makefile) that using the Lapack routine dgesv solves the following system of equations:
\begin{eqnarray}
3 x_{1} + 2 x_{2} & = & 2 \\
2 x_{1} + 6 x_{2} & = & -8
\end{eqnarray}
Documentation for the library can be found [[here|http://www.netlib.org/lapack/double/dgesv.f]]. This document assumes FORTRAN language. If you use C, call dgesv_(&n,&nrhs,a,&lda,ipiv,b,&ldb,&info) (note the underscore). The arguments and usage information are the same in both languages.

!!! 2. Scaling analysis of LAPACK routine dgesv
We wish to examine how the execution time grows with $n$, the order of the system of equations to solve. In order to do so,
# Create a large matrix $A$ with random entries. CALL srand(int seed) (seed is an arbitrary integer) will initialize the random sequence, and rand() returns real random numbers uniformly distributed between 0 and 1. Create a matrix of size $(n,n)$ with random entries, and also a right hand side vector $b$. To make the problem better conditioned, add to the diagonal elemens of $A$ some positive quantity (e.g. 1).
# Solve the linear system of equations for increasing values of $n$. Use the time command to obtain the execution time. For example, if a.out is the executable, type instead 
{{{
time a.out
}}}
This will give you the execution time. Plot in a log log scale the execution time with $n$ and examine how the time scales with $n$ for large $n$.
Although simple, this method will also include the time involved in the other tasks in your code. For example, it will include the time spent setting up the matrix A. A more accurate timing can be obtained by using the function cpu_time (in fortran):
{{{
CALL cpu_time(t0)
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE(6,*) 'Elapsed time (in seconds): ',t1-t0
}}}
In C you would [[use|https://www.gnu.org/software/libc/manual/html_node/CPU-Time.html]]:
{{{
#include <time.h>
clock_t start, end;
double cpu_time_used;

start = clock();
… /* Do the work. */
end = clock();
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
}}}


!!! 3. Zero of function with spline interpolation
We have found experimentally the pressure $p$ as a function of the molar volume $v$ of a binary mixture at constant temperature and composition (in some arbitrary units). The results are shown in the following table.
|v | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 |
| p | 4.7 | 3.5 | 3.0 | 2.7 | 3.0 | 2.4 |
# Write your own code to find a cubic spline interpolation polynomial $p=p(v)$ of this data in the interval given. Plot the data along with the spline interpolation on the same graph.
# Use the spline interpolation to find the value of the volume at which $p=3.25$. Do it approximately on the graph first.
# Compute the value of the volume at which $p=3.25$ by solving the equation $f(v)=p(v)-3.25=0$ by using the ~Newton-Raphson method ''WITH'' the interpolated spline function.
# Additional information
## It is not necessary for the solution of this problem, but there is [[documentation|http://mathworld.wolfram.com/CubicSpline.html]] about the definition and calculation of spline interpolating polynomials.
## If you use FORTRAN, you can use the routines spline.f and splint.f provided in the class github repository.
## If you use C, you can use the routines spline.c and splint.c and header files nrutil.h, nrutils.c provided in the class github repository.
{{justifyright{
''Due Tuesday, September 25, 2018''
}}}

!!! 1. Bessel function approximation
Bessel functions are transcendental functions that can be computed by summing a series. However, there are recursion formulae that make the computation more efficient, although, as seen in class, recursion can introduce its own numerical problems. Use the following recursion formulae to compute the Bessel function $J_{l}$ with $l = 0, 1, \ldots, 8$ and plot the results for $x \in (0,12)$. The recursion formulae ("upwards" or "downwards") are \begin{eqnarray}
J_{l+1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l-1}(x) \\
J_{l-1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l+1}(x) \\
{\rm Note:} \;\; J_{0}(x) = \frac{\sin x}{x} & \; & J_{1}(x) = \frac{\sin x - x \cos x}{x^{2}}
\end{eqnarray}  The following table lists accurate values of selected Bessel functions, and can be used to help evaluate the accuracy of your approximations,
|! x |! $J_{3}$ |! $J_{5}$ |! $J_{8}$ |
| 0.10 | 9. 518519719 10$^{-6}$ | 9.616310231 10$^{-10}$ | 2.901200102 10$^{-16}$ |
| 1.00 | 9.006581118 10$^{-3}$ | 9.256115862 10$^{-5}$ | 2.826498802 10$^{-8}$ |
| 10.0 | -3.949584498 10$^{-2}$ | -5.553451162 10$^{-1}$ | 1.255780236 10$^{-1}$ |
# Write a code that computes the Bessel functions required by using the upward recursion.
# Compare your solution to the table. Argue how errors might have propagated and what effect this may have on the accuracy of the solution.
# Chose arbitrary starting values for $J_{L}(x)$ and $J_{L-1}(x)$ with large $L$ and use the iteration downward. Then use the ratio between the analytical and recursive values for $J_{0}(x)$ to normalize the value for the desired $J_{l } (x)$. Compare the results to the upward iteration for the Bessel functions asked above.

!!! 2. Solution of simple linear system of equations (LAPACK)
Write a code (and makefile) that using the Lapack routine dgesv solves the following system of equations:
\begin{eqnarray}
3 x_{1} + 2 x_{2} & = & 2 \\
2 x_{1} + 6 x_{2} & = & -8
\end{eqnarray}
Documentation for the library can be found [[here|http://www.netlib.org/lapack/double/dgesv.f]]. This document assumes FORTRAN language. If you use C/C++, call dgesv_(&n,&nrhs,a,&lda,ipiv,b,&ldb,&info) where the length of ''array'' a is $lda*n$ (note also the underscore). Documentation [[here|http://www.physics.utah.edu/~detar/phys6720/handouts/lapack.html]]. The arguments and usage information are the same in both languages. In C++ keep in mind the extra statement needed to indicate that dgesv has been compiled in C.

!!! 3. Scaling analysis of LAPACK routine dgesv
We wish to examine how the execution time grows with $n$, the order of the system of equations to solve. In order to do so, write a program that
# Creates a large real matrix of order $n$, $A(n,n)$, with random entries. In FORTRAN, use CALL srand(seed) (where seed is an arbitrary integer) to initialize the random sequence. In C, use the function void srand(unsigned int seed) instead. Then, in both languages, the function rand() returns real random numbers uniformly distributed between 0 and 1. Create a matrix of size $(n,n)$ with random entries, and also a right hand side vector $b$. To make the problem better conditioned, add to the diagonal elements of $A$ some positive quantity (e.g. 1).
# Solve the linear system of equations for increasing values of $n$. Use the time command to obtain the execution time. For example, if a.out is the executable, type instead 
{{{
time a.out
}}}
This will give you the execution time. Plot in a log log scale the execution time with $n$ and examine how the time scales with $n$ for large $n$.
Although simple, this method will also include the time involved in the other tasks in your code. For example, it will include the time spent setting up the matrix A. A more accurate timing can be obtained by using the function cpu_time (in fortran):
{{{
CALL cpu_time(t0)
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE(6,*) 'Elapsed time (in seconds): ',t1-t0
}}}
In C you would [[use|https://www.gnu.org/software/libc/manual/html_node/CPU-Time.html]]:
{{{
#include <time.h>
clock_t start, end;
double cpu_time_used;

start = clock();
… /* Do the work. */
end = clock();
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
}}}
{{justifyright{
''Due Tuesday, October 8, 2014''
}}}
!!!!1. Solution of Poisson equation
Poisson's equation relates the electrostatic potential $V(x)$ to a charge distribution $\rho(x)$. In arbitrary units, and in two spatial dimensions we write,
\begin{equation}
\frac{\partial^{2} V(x,y)}{\partial x^{2}} + \frac{\partial^{2} V(x,y)}{\partial y^{2}} = - \rho(x,y)
\end{equation}

By a specific discretization scheme, the equation is written as ($\Delta x = \Delta y$):
\begin{equation}
\frac{1}{\Delta x^{2}} \left( V_{i+1,j} - 2 V_{i,j} + V_{i-1,j} + V_{i,j+1} - 2 V_{i,j} + V_{i,j-1} \right) = - \rho_{i,j} \quad i = 1, ..., N; ~ j = 1, ... , M
\end{equation}
In this form, given a desity $\rho$, and appropriate boundary conditions (e.g., values of $V$ on the edges of the rectangular region implict in the discretization given), solving the Poisson equation is recast as a linear problem: $Ax=b$. $A$ is a sparse matrix, $x$ is a vector of length $M \times N$, and $b$ contains the given density $\rho$.

Write a code that would use an iterative method of your choice to solve Poisson's equation. Imagine that the charge density is given, and that the boundary consitions are also specified.

!!!!2. Numerical solution of Shrodinger's equation
Consider a one dimensional potential well that is a truncated harmonic oscillator potential:
\begin{equation}
V(x) = \frac{1}{2} x^{2} \quad 0 < x < 1, \quad\quad  V = \infty \quad x=0,1
\end{equation}
Take as basis functions,
\begin{equation}
u_{i}(x) = \sin (q_{i} x) \quad q_{i} = i \pi
\end{equation}

Consider the Hamiltonian (in arbitrary units),
\begin{equation}
H = - \frac{d^{2} }{dx^{2}} + V(x) 
\end{equation}

We wish to numerically find the energy levels of this system. First, compute the necessary matrix elements:
\begin{equation}
H_{ji} = \langle u_{j} | H | u_{i} \rangle = \int_{0}^{1} u_{j}(x) H u_{i}(x) dx 
\end{equation}
Given the normalization,
\begin{equation}
S_{ji} = \langle u_{j} | u_{i} \rangle = \int_{0}^{1} u_{j}(x) u_{i}(x) dx
\end{equation}  
then, Schrodinger's equation turns into an eigenvalue problem of the form:
\begin{equation}
H \alpha = E S \alpha
\end{equation}
where $\alpha$ are the components of $\psi$ in the basis given $\psi = \sum_{i} \alpha_{i} u_{i}(x)$.

Use a Lapack routine to find the energy levels of this potential. Make a judicious choice about the number of eigenfunctions you use in the diagonalization.
{{justifyright{
''Due Friday, October 7, 2016''
}}}
!!! 1. The Sucessive Over Relaxation method (SOR)
# Show that the ~Gauss-Seidel iterative method seen in class can be written as $$ x_{i}^{(k+1)} = \frac{- \sum_{j=1}^{i-1} A_{ij} x_{j}^{(k+1)}-\sum_{j=i+1}^{n} A_{ij} x_{j}^{(k)} + b_{i}}{A_{ii}}, \quad i=1, \ldots, n $$
where the residual to be made as small as possible has been defined as $r_{i}^{(k)} = x_{i}^{(k+1)} - x_{i}^{(k)}$.
# "Over relaxation" is introduced in order to speed up convergence of the iteration, and consists in introducing a new iteration that overcorrects the residual,$$ x_{i}^{(k+1)} = x_{i}^{(k)} + \omega r_{i}^{(k)}$$ where $\omega$ is known as the relaxation parameter. Show that if we write, in matrix form, $x^{(k+1)} = B x^{(k)} + C$, then for this new method we have for the iteration  $$ B_{\rm SOR} = \left( \mathbb{1} + \omega D^{-1}L \right)^{-1} \left[ (1-\omega)\mathbb{1} - \omega D^{-1}U \right]$$
with the same notation used in class for $D, L$ and $U$.
# Show that the determinant of $B$ is $\det (B) = (1-\omega)^{n}$.
# To estimate the spectral radius (the region of convergence of the method), assume that $B$ can be diagonalized, with eigenvalues $\lambda_{1}, \lambda_{2}, \ldots \lambda_{n}$. With this representation show that $\| B \| \ge |1-\omega|$. What is the range of $\omega$ over which this procedure will converge ?
!!! 2. Solution of a system of equations by iterative methods
Consider the system of equations $A \mathbf{x} = \mathbf{b}$ where $$ A = \left( \begin{array}{cccc} 4 & -1 & -1 & 0 \\ -1 & 4 & 0 & -1 \\ -1 & 0 & 4 & -1 \\ 0 & -1 & -1 & 4 \end{array} \right), \quad\quad b = \left( \begin{array}{c} 1 \\ 2 \\ 0 \\1 \end{array} \right). $$
Write a code that solves this system by using (i) the Jacobi method, (ii) the ~Gauss-Seidel method, and (iii) the SOR method. From the output of the iterations, compare empirically the rate of convergence of the three methods.
!!!3. Solution of Poisson equation in two dimensions
Poisson's equation relates the electrostatic potential $V(x,y)$ to a charge distribution $\rho(x,y)$. In arbitrary units, and in two spatial dimensions we write,
\begin{equation}
\frac{\partial^{2} V(x,y)}{\partial x^{2}} + \frac{\partial^{2} V(x,y)}{\partial y^{2}} = - \rho(x,y)
\end{equation}

By some specific discretization scheme, the equation is written as (with $V_{ij} = V(x_{i},y_{j})$, $\Delta x = \Delta y$):
\begin{equation}
\frac{1}{\Delta x^{2}} \left( V_{i+1,j} - 2 V_{i,j} + V_{i-1,j} + V_{i,j+1} - 2 V_{i,j} + V_{i,j-1} \right) = - \rho_{i,j} \quad i = 1, ..., N; ~ j = 1, ... , M
\end{equation}
In this form, given a density $\rho$, and appropriate boundary conditions (e.g., values of $V$ on the edges of the rectangular region implicit in the discretization given), you can think of solving the Poisson equation as solving a linear problem: $Ax=b$. $A$ is a sparse matrix, $x$ is a vector of length $M \times N$, and $b$ contains the given density $\rho$. Alternatively, you may also directly write an iterative procedure to solve the discrete system of $N \times M$ equations for the potential $V_{ij}$.

Write a code that would use an iterative method of your choice to solve Poisson's equation. Imagine that the charge density is given, and that the boundary conditions are also specified.
{{justifyright{
''Due Tuesday, October 2, 2018''
}}}

!!! Question 1. Spline interpolation with explicit solution of a tridiagonal system
Spline polynomials are smooth, piecewise polynomials that are extensively used in Physics to interpolate a sample of discrete points.

We have found experimentally the pressure $p$ as a function of the molar volume $v$ of a binary mixture at constant temperature and composition (in some arbitrary units). The results are shown in the following table.
|v | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 |
| p | 4.7 | 3.5 | 3.0 | 2.7 | 3.0 | 2.4 |
# Write your own code to find a cubic spline interpolation polynomial $p=p(v)$ of this data in the interval given. Plot the data along with the spline interpolation on the same graph.
# Use the spline interpolation to find the value of the volume at which $p=3.25$. Do it approximately on the graph first.
# Compute the value of the volume at which $p=3.25$ by solving the equation $f(v)=p(v)-3.25=0$ by using the ~Newton-Raphson method ''WITH'' the interpolated spline function.
# Additional information
## Documentation about spline polynomials can be found at Wolfram [[MathWorld|http://mathworld.wolfram.com/CubicSpline.html]]. Eq. (1) defines the collection of polynomials, and Eqs. (6-9) define the mapping between these coefficients and the set $D_{i}$. The latter is the one that is actually calculated in practice by solving the tridiagonal system in Eq. (18).
## Use the algorithm given in class to solve tridiagonal equations. 

!!! Question 2. Spline interpolation. Routines provided
Rewrite the code of question 1 but now used spline interpolation routines used by others.
# If you use FORTRAN, you can use the routines spline.f to compute the spline coefficients and seval.f to evaluate the polynomials at interpolated points. The routines are provided in the class github repository.
# If you use C++, you can use the functions spline.c to compute the polynomial coefficients, splint.c to evaluate the polynomials obtained. You will also need to add in the compilation nrutil.c and the header file nrutil.h. All of them are provided in the class github repository.

!!! Question 3. Numerical solution of Schroedinger's equation
Consider a one dimensional potential well that is a truncated harmonic oscillator potential:
\begin{equation}
V(x) = \frac{1}{2} x^{2} \quad 0 < x < 1, \quad\quad  V = \infty \quad x=0,1
\end{equation}
Take as basis functions,
\begin{equation}
u_{i}(x) = \frac{1}{\sqrt{2}} \sin (q_{i} x) \quad q_{i} = i \pi
\end{equation}

Consider the Hamiltonian (in arbitrary units),
\begin{equation}
H = - \frac{d^{2} }{dx^{2}} + V(x) 
\end{equation}

We wish to numerically find the energy levels of this system. First, compute the necessary matrix elements by using the trapezoidal rule,
\begin{equation}
H_{ji} = \langle u_{j} | H | u_{i} \rangle = \int_{0}^{1} u_{j}(x) H u_{i}(x) dx .
\end{equation}
Schroedinger's equation turns into an eigenvalue problem of the form:
\begin{equation}
H \alpha = E \alpha
\end{equation}
where $\alpha$ are the components of $\psi$ in the basis given $\psi = \sum_{i} \alpha_{i} u_{i}(x)$.

Use a Lapack routine to find the energy levels of this potential. Make a judicious choice about the number of eigenfunctions you use in the diagonalization.
{{justifyright{
''Due Tuesday, October 14, 2014''
}}}
# Write code that verifies the relatively faster convergence of a numerical derivative computed by central differences compared to a derivative computed by forward differences. Evaluate the lowest order forward and central difference approximation to the derivative of $f(x) = \tan (x)$ in the interval $[-2,2]$ for $n$ equally spaced points. Compare the results to the exact values of the derivative. For a given $n$ find the maximum error on the interval from each method, defining the error as the absolute value of the ratio of the difference between numerical and analytic derivatives. Then make a log-log plot of this error versus $n$ over at least one order of magnitude in $n$ (say from $n = 100$ to $n=5000$) to show that the forward difference estimate converges linearly with $n$, while the central difference estimate converges as $n^{2}$.
# Numerical intergation over an unbounded domain. The integration methods discussed in class consider definite integrals over a finite domain $(a,b)$. However, it is possible to either transform the integrands or to analyze convergence with the limits of integration to evaluate integrals over unbounded domains. Propose a suitable plan, and write the associate code to evaluate the following integrals for $s = 0.5, 0.6, 0.7, ..., 3.0$. 
## $$\int_{0}^{\infty} \left( x^{3} + sx \right)^{-1/2} dx,$$  
## $$\int_{0}^{\infty} \left( x^{2} + 1 \right) ^{-1/2} e^{-sx} dx$$ 
# Numerical integration of a predator-prey model. Consider two species of fish that coexist in the same lake. Let $x$ be the number of small fish of type A which only eat grass from the bottom, and $y$ the number of large fish of type B which only eat small fish of type A, but no grass. Alone, the population of A might be expected to grow exponentially, whereas B alone (no A to eat) would die exponentially. However, the two species interact: whenever a fish B encounters a fish A, it will eat it. This interaction contributes to the increase in the population of B and a decrease in A, both rates being proportional to the quantity $xy$. A general mathematical model of this system is $$ \frac{dx}{dt} = \epsilon_{A} x - \gamma_{A} xy $$
$$ \frac{dy}{dt} = - \epsilon_{B} y + \gamma_{B} xy $$  where $\epsilon_{A}$ is the growth rate of A alone, $\epsilon_{B}$ the extinction of B alone, and $\gamma_{A}$ and $\gamma_{B}$ the interaction coefficients. All four are assumed to be constant. Write a code that integrates this system of differential equations with $\epsilon_{A} = \epsilon_{B} = \gamma_{A} = \gamma_{B} =1$, and initial conditions at $t=0$ $x=0.1$ and $y=0.1$. 
## Integrate the system up to $t=10$.
## Plot $y$ versus $x$ and provide an interpretation of the resulting graph.

[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/PHYS4041_Homework4.pdf]].
{{justifyright{
''Due Friday, October 14, 2016''
}}}
!!!1. Numerical solution of Shrodinger's equation
Consider a one dimensional potential well that is a truncated harmonic oscillator potential:
\begin{equation}
V(x) = \frac{1}{2} x^{2} \quad 0 < x < 1, \quad\quad  V = \infty \quad x=0,1
\end{equation}
Take as basis functions,
\begin{equation}
u_{i}(x) = \sin (q_{i} x) \quad q_{i} = i \pi
\end{equation}

Consider the Hamiltonian (in arbitrary units),
\begin{equation}
H = - \frac{d^{2} }{dx^{2}} + V(x) 
\end{equation}

We wish to numerically find the energy levels of this system. First, compute the necessary matrix elements:
\begin{equation}
H_{ji} = \langle u_{j} | H | u_{i} \rangle = \int_{0}^{1} u_{j}(x) H u_{i}(x) dx 
\end{equation}
Given the normalization,
\begin{equation}
S_{ji} = \langle u_{j} | u_{i} \rangle = \int_{0}^{1} u_{j}(x) u_{i}(x) dx
\end{equation}  
then, Schrodinger's equation turns into an eigenvalue problem of the form:
\begin{equation}
H \alpha = E S \alpha
\end{equation}
where $\alpha$ are the components of $\psi$ in the basis given $\psi = \sum_{i} \alpha_{i} u_{i}(x)$.

Use a Lapack routine to find the energy levels of this potential. Make a judicious choice about the number of eigenfunctions you use in the diagonalization.

!!!2. Convergence of numerical differentiation
Write code that verifies the relatively faster convergence of a numerical derivative computed by central differences compared to a derivative computed by forward differences. Evaluate the lowest order forward and central difference approximation to the derivative of $f(x) = \tan (x)$ in the interval $[-2,2]$ for $n$ equally spaced points. Compare the results to the exact values of the derivative. For a given $n$ find the maximum error on the interval from each method, defining the error as the absolute value of the ratio of the difference between numerical and analytic derivatives. Then make a log-log plot of this error versus $n$ over at least one order of magnitude in $n$ (say from $n = 100$ to $n=5000$) to show that the forward difference estimate converges linearly with $n$, while the central difference estimate converges as $n^{2}$.

!!!3. Numerical integration over an unbounded domain. 
The integration methods discussed in class consider definite integrals over a finite domain $(a,b)$. However, it is possible to either transform the integrands or to analyze convergence with the limits of integration to evaluate integrals over unbounded domains. Propose a suitable plan, and write the associate code to evaluate the following integrals for $s = 0.5, 0.6, 0.7, ..., 3.0$. 
## $$\int_{0}^{\infty} \left( x^{3} + sx \right)^{-1/2} dx,$$  
## $$\int_{0}^{\infty} \left( x^{2} + 1 \right) ^{-1/2} e^{-sx} dx$$ 
{{justifyright{
''Due Tuesday, October 9, 2018''
}}}

!!! 1. Solution of a system of equations by iterative methods
Consider the system of equations $A \mathbf{x} = \mathbf{b}$ where $$ A = \left( \begin{array}{cccc} 4 & -1 & -1 & 0 \\ -1 & 4 & 0 & -1 \\ -1 & 0 & 4 & -1 \\ 0 & -1 & -1 & 4 \end{array} \right), \quad\quad b = \left( \begin{array}{c} 1 \\ 2 \\ 0 \\1 \end{array} \right). $$
Write a code that solves this system by using (i) the Jacobi method, and (ii) the ~Gauss-Seidel method. From the output of the iterations, compare empirically the rate of convergence of the two methods.

!!! 2. Solution of an integral equation
We need to find the solution of the following integral equation for $f(x)$ in the interval $(0,1)$,$$
\int_{0}^{1} dx e^{-(x-y)^{2}/\epsilon^{2}} \frac{df(x)}{dx} = 2 \epsilon \sqrt{\pi} ( y - 3 y^{2} + 2 y^{3}),
$$
where the function $f(x)$ needs to satisfy the boundary conditions $f(0)=f(1)=0$, and $\epsilon = 0.01$.

Discretize $x_{j}= j h, j = 0, \ldots n$ (take $n = 101$) for some $h$, so that $x_{0} = 0$ and $x_{n} = 1$, and $y_{i} = i h, i = 0, \ldots n$. Consistent with this discretization, one defines $f_{j} = f(x_{j})$. The boundary conditions are then $f_{0} = f_{n} = 0$.

# Write the discretized version of the integral equation by using the trapezoidal rule for the integral and central differences for the derivative at the inner points, the forward difference at $x=0$ and the backward difference at $x=1$.
# Use the discretized version to write the code that will find the solution $f(x)$ in $(0,1)$.
{{justifyright{
''Due Tuesday, November 4/11, 2014''
}}}
!!!!1. ~Lax-Wendroff method. Due November 4.
Write a program that uses the ~Lax-Wendroff method to solve the PDE:
$$ \frac{\partial u}{\partial t} = - v \frac{\partial u}{\partial x}$$
with $v$ a given constant. As given in class, the method as applied to this equation is
$$u_{j}^{n+1}=u_{j}^{n}-\frac{v \Delta t}{\Delta x} \left( u_{j+1/2}^{n+1/2} - u_{j-1/2}^{n+1/2} \right)$$
with
$$u_{j+1/2}^{n+1/2} = \frac{1}{2} \left( u_{j+1}^{n} + u_{j}^{n} \right) - \frac{v \Delta t}{2 \Delta x} \left( u_{j+1}^{n} - u_{j}^{n} \right) $$
and a similar definition for $u_{j-1/2}^{n+1/2}$. Space is discretized by $\Delta x$ and $j$ is the spatial index, whereas time has been discretized by $\Delta t$ with corresponding index $n$.

Let $a = |v|\Delta t/\Delta x \le 1$ denote the Courant condition for stability. 

Consider as initial condition $u(x,t=0) = \sin (2 \pi x)$ for $0 \le x \le 1$. If your discretization is such that $x_{0} = 0$ and $x_{N} = 1$, have the code solve the equation with "periodic boundary conditions": $u_{N} = u_{0}$ at each time (and should you need it $u_{N+1} = u_{1}$, etc.). For the purposes of the homework, choose $N = 256$ and $v = 1$.

Run your code with $a = 0.99$. Since $v = 1$, the solution  is advected in such a way that the solution should exactly lie on top of the initial condition for times $t = 1, 2, 3, ...$. Confirm this graphically for $t = 1, 2$.

Run the code until $t = 0.1$ using both $a = 0.99$ and $a = 1.1$. Plot the solutions and comment on any differences.

!!!!2. Steady state temperature field in a porous medium. Due November 11.
The temperature distribution in a porous rock in the presence of ground water flow directed along the $x$ axis is described by the equation
$$ v \frac{\partial T}{\partial x} = \kappa \nabla^{2} T$$
where $v$ is the constant water velocity and $\kappa$ is the thermal diffusivity coefficient. The domain of solution is the rectangle $0 < xv/\kappa < 30$ $0 < yv/\kappa < 10$.

Compute the temperature distribution using the Alternating Direction Implicit method combined with a finite difference method of your choice. The boundary conditions $T = 65$ at $x=0$, $T = 25$ at $x = 30 \kappa/v$, $T = 25$ at $y = 10 \kappa/v$ and $\partial T/\partial y = 0$ at $y=0$.

!!!! Notes on the ADI method
It is useful to be able to write the system of equations of Problem 2 as the usual system $Ax = b$ where both $x$ and $b$ are arrays. Consider for example the two dimensional Poisson equation
$$ T_{i+1 j} - 2 T_{ij} + T_{i11 j} + T_{i j+1} - 2 T_{ij} + T_{i j-1} = - \Delta x^{2} \rho_{ij}$$
for $i = 1, ... , N$ and $j = 1, ..., M$.

We can now construct a long vector $\mathbf{v}$ of length $N \times M$ by linking together the @@color(red):rows@@ of the matrix $T_{ij}$:
$$ v_{r} = T_{ij} \quad \quad r = (i-1)M + j$$
This mapping can be inverted as
$$ i = 1+ {\rm INT}\left( \frac{r-1}{M} \right) \quad\quad j = 1 + (r-1){\rm MOD}(M)$$ 
The Poisson equation then transforms to
$$ v_{r-M} + v_{r-1} - 4 v_{r} + v_{r+1} + v_{r+M} = - \Delta x^{2} \rho_{r}$$ which can now written in the standard form $Ax = b$.

For the ADI method we need to construct another long vector $\mathbf{w}$ by linking together the @@color(red):columns@@ of the matrix $T_{ij}$:
$$ w_{s} = T_{ij} \quad\quad s = (j-1)N + i $$
which can in turn be inverted as
$$ j = 1 + {\rm INT}\left( \frac{s-1}{N} \right) \quad\quad i = 1 + (s-1){\rm MOD}(N) $$
With the aid of @@color(red):both@@ vectors, the Poisson equation can be written as
$$ w_{s+1} - 2w_{s} + w_{s-1} + v_{r+1} - 2 v_{r} + v_{r-1} = - \Delta x^{2} \rho_{ij}$$
which one can write in matrix form as
$$ A_{1} \mathbf{v} + A_{2} \mathbf{w} = b$$
so that $A_{1}$ only acts on the rows of $T_{ij}$ and $A_{2}$ on the columns. We also define the reordering matrix $U$ as
$$ \mathbf{w} = U \mathbf{u}$$

In this notation, the ADI iteration can be written as
\begin{eqnarray}
\left( A_{1} + I \right) \mathbf{v}^{n+1/2} & = & b - (A_{2} \mathbf{w}^{n} - \mathbf{v}^{n} ) \\
\mathbf{w}^{n+1/2} & = & U \mathbf{v}^{n+1/2} \\
\left( A_{2} + I \right) \mathbf{w}^{n+1} & = & b - (A_{1} \mathbf{v}^{n+1/2} - \mathbf{w}^{n+1/2} )
\end{eqnarray}
where $I$ is the identity matrix.
{{justifyright{
''Due Friday, November 4, 2016''
}}}

!!! ~Lax-Wendroff method.
Write a program that uses the ~Lax-Wendroff method to solve the PDE:
$$ \frac{\partial u}{\partial t} = - v \frac{\partial u}{\partial x}$$
with $v$ a given constant. As given in class, the method as applied to this equation is
$$u_{j}^{n+1}=u_{j}^{n}-\frac{v \Delta t}{\Delta x} \left( u_{j+1/2}^{n+1/2} - u_{j-1/2}^{n+1/2} \right)$$
with
$$u_{j+1/2}^{n+1/2} = \frac{1}{2} \left( u_{j+1}^{n} + u_{j}^{n} \right) - \frac{v \Delta t}{2 \Delta x} \left( u_{j+1}^{n} - u_{j}^{n} \right) $$
and a similar definition for $u_{j-1/2}^{n+1/2}$. Space is discretized by $\Delta x$ and $j$ is the spatial index, whereas time has been discretized by $\Delta t$ with corresponding index $n$.

Let $a = |v|\Delta t/\Delta x \le 1$ denote the Courant condition for stability. 

Consider as initial condition $u(x,t=0) = \sin (2 \pi x)$ for $0 \le x \le 1$. If your discretization is such that $x_{0} = 0$ and $x_{N} = 1$, have the code solve the equation with "periodic boundary conditions": $u_{N} = u_{0}$ at each time (and should you need it $u_{N+1} = u_{1}$, etc.). For the purposes of the homework, choose $N = 256$ and $v = 1$.

Run your code with $a = 0.99$. Since $v = 1$, the solution  is advected in such a way that the solution should exactly lie on top of the initial condition for times $t = 1, 2, 3, ...$. Confirm this graphically for $t = 1, 2$.

Run the code until $t = 0.1$ using both $a = 0.99$ and $a = 1.1$. Plot the solutions and comment on any differences.
{{justifyright{
''Due Tuesday, October 30, 2018''
}}}
!!! 1. Eigenvalue problem for a differential equation
Write a code that computes the smallest eigenvalue $\lambda$ for which the problem $$
\frac{d}{dx} \left[ (1+x^{2}) \frac{dy}{dx} \right] + \lambda y = 0, \quad\quad y(-1)=y(1)=0
$$
has a nontrivial solution.

In order to do it, you can use one of the following two methods. Either follow the method described in class in which a new variable $z(x)=\lambda$ is introduced, with $dz/dx = 0$, or user a matrix method in the following way: discretize the equation as, $$ 
\frac{d}{dx} \left[ p(x) \frac{dy}{dx} \right] = \frac{1}{h} \left[ p_{n+1/2}\left( \frac{y_{n+1}-y_{n}}{h} \right) - p_{n-1/2} \left( \frac{y_{n}-y_{n-1}}{h} \right) \right] + c_{1}(x) h^{2} + c_{2} h^{4} + \ldots ,
$$
where we have also included the asymptotic error of the discretization. Once the differential equation is discretized, you can find the smallest eigenvalue of the resulting matrix.



!!! 2. ~Lax-Wendroff method
Write a program that uses the ~Lax-Wendroff method to solve the PDE:
$$ \frac{\partial u}{\partial t} = - v \frac{\partial u}{\partial x}$$
with $v$ a given constant. As given in class, the method as applied to this equation is
$$u_{j}^{n+1}=u_{j}^{n}-\frac{v \Delta t}{\Delta x} \left( u_{j+1/2}^{n+1/2} - u_{j-1/2}^{n+1/2} \right)$$
with
$$u_{j+1/2}^{n+1/2} = \frac{1}{2} \left( u_{j+1}^{n} + u_{j}^{n} \right) - \frac{v \Delta t}{2 \Delta x} \left( u_{j+1}^{n} - u_{j}^{n} \right) $$
and a similar definition for $u_{j-1/2}^{n+1/2}$. Space is discretized by $\Delta x$ and $j$ is the spatial index, whereas time has been discretized by $\Delta t$ with corresponding index $n$.

Let $a = |v|\Delta t/\Delta x \le 1$ denote the Courant condition for stability. 

Consider as initial condition $u(x,t=0) = \sin (2 \pi x)$ for $0 \le x \le 1$. If your discretization is such that $x_{0} = 0$ and $x_{N} = 1$, have the code solve the equation with "periodic boundary conditions": $u_{N} = u_{0}$ at each time (and should you need it $u_{N+1} = u_{1}$, etc.). For the purposes of the homework, choose $N = 256$ and $v = 1$.

Run your code with $a = 0.99$. Since $v = 1$, the solution  is advected in such a way that the solution should exactly lie on top of the initial condition for times $t = 1, 2, 3, ...$. Confirm this graphically for $t = 1, 2$.

Run the code until $t = 0.1$ using both $a = 0.99$ and $a = 1.1$. Plot the solutions and comment on any differences.

{{justifyright{
''Due Friday, November 11, 2016''
}}}

!!! Steady state temperature field in a porous medium.
The temperature distribution in a porous rock in the presence of ground water flow directed along the $x$ axis is described by the equation
$$ v \frac{\partial T}{\partial x} = \kappa \nabla^{2} T$$
where $v$ is the constant water velocity and $\kappa$ is the thermal diffusivity coefficient. The domain of solution is the rectangle $0 < xv/\kappa < 30$ $0 < yv/\kappa < 10$.

In order to find the solution, we will solve the time dependent convection-diffusion equation instead $$
\frac{\partial T}{\partial t} = - v \frac{\partial T}{\partial x} + \kappa \left( \frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} \right), 
$$
and iterate it in time until it reaches a steady state $\partial T/\partial t = 0$, at which point we have solved the required equation.

By using forward differences for the time derivative, central differences for the spatial derivatives, and the Alternatig Direction Implicit method to consider a semi-implicit solution in the two spatial directions, obtain the temperature distribution given the following boundary conditions are $T = 65 \;$ at  $\; x=0$, $T = 25 \;$  at $\; x = 30 \kappa/v$, $T = 25 \;$ at $\; y = 10 \kappa/v \;$ and $\; \partial T/\partial y = 0 \;$ at $\; y=0$. Note, you can use $\kappa = v = 1$.
{{justifyright{
''Due Tuesday, November 6, 2018''
}}}
!!! Steady state temperature field in a porous medium.
The temperature distribution in a porous rock in the presence of ground water flow directed along the $x$ axis is described by the equation
$$ v \frac{\partial T}{\partial x} = \kappa \nabla^{2} T$$
where $v$ is the constant water velocity and $\kappa$ is the thermal diffusivity coefficient. The domain of solution is the rectangle $0 < xv/\kappa < 30$ $0 < yv/\kappa < 10$.

In order to find the solution, we will solve the time dependent convection-diffusion equation instead $$
\frac{\partial T}{\partial t} = - v \frac{\partial T}{\partial x} + \kappa \left( \frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} \right), 
$$
and iterate it in time until it reaches a steady state $\partial T/\partial t = 0$, at which point we have solved the required equation.

By using forward differences for the time derivative, central differences for the spatial derivatives, and the Alternatig Direction Implicit method to consider a semi-implicit solution in the two spatial directions, obtain the temperature distribution given the following boundary conditions are $T = 65 \;$ at  $\; x=0$, $T = 25 \;$  at $\; x = 30 \kappa/v$, $T = 25 \;$ at $\; y = 10 \kappa/v \;$ and $\; \partial T/\partial y = 0 \;$ at $\; y=0$. Note, you can use $\kappa = v = 1$.
!!!! Notes on the ADI method
As discussed in class, it is useful to be able to write the system of equations as the usual system $Ax = b$ where both $x$ and $b$ are arrays. We discretize the two dimensional domains as follows $i = 1, ... , N$ and $j = 1, ..., M$. In the first half step (in the same sequence of two steps as discussed in class) we construct a long vector $\mathbf{v}$ of length $N \times M$ by linking together the @@color(red):columns@@ of the matrix $T_{ij}$:
$$ v_{l} = T_{ij} \quad \quad l = (j-1)N + i$$
This mapping can be inverted as
$$ i = 1+ {\rm INT}\left( \frac{l-1}{N} \right) \quad\quad j = 1 + (l-1){\rm MOD}(N)$$ 

For the second half step, we need to construct another long vector $\mathbf{w}$ by linking together the @@color(red):rows@@ of the matrix $T_{ij}$:
$$ w_{m} = T_{ij} \quad\quad m = (i-1)M + j$$
which can in turn be inverted as
$$ j = 1 + {\rm INT}\left( \frac{m-1}{M} \right) \quad\quad i = 1 + (m-1){\rm MOD}(M) $$.
{{justifyright{
''Due Tuesday, November 25, 2014''
}}}
!!Molecular Dynamics simulation
Write a Molecular Dynamics simulation code with the following features:
* Consider a fluid system in a ''two dimensional'' square box of size $L \times L$.
* Particles interact via a standard ~Lennard-Jones potential $$V(r) = 4 \epsilon \left( \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r}\right)^{6} \right)$$
where $\epsilon$ and $\sigma$ are constants and assumed known.
* Use the velocity-Verlet algorithm for time integration.
* The system will be held at constant temperature $T$. In order to accomplish that, special boundary conditions will be assumed at the walls of the box (a thermostat at constant temperature). When a particle collides with any of the walls, the tangential component of the velocity of the particle is conserved, and the incident normal velocity $v_{n}$ is updated into an outgoing normal velocity $v_{n}^{\prime}$ given by [ [[van Beijeren|http://arxiv.org/abs/1411.2983]], 2014]: $$v_{n}^{\prime} = \sqrt{ - \frac{2}{\beta m} \ln \left( 1 - {\rm exp} \left( - \frac{\beta m v_{n}^{2}}{2} \right) \right) }$$ Or course, $\beta = 1 /k_{B}T$, and $m$ is the mass of the particles, all assumed identical.
* Make sure you use appropriate dimensionless variables.

Run the code and compute the internal energy, pressure and pair correlation function at some reasonable temperature.

!!!Notes on dimensionless units
Choose $\sigma$ as the unit of length: $r^{\prime}=r/\sigma$ is the dimensionless distance, and $L^{\prime}=L/\sigma$ the dimensionless lateral dimension of the box (e.g. 10). If we introduce a dimensionless potential energy $V^{\prime}=V/\epsilon$ and a dimensionless time $t^{\prime}=t/\tau$, then Newton's Law transforms to
$$ \frac{m \sigma^{2}}{\tau^{2} \epsilon} \frac{d^{2} r^{\prime}}{dt^{\prime \; 2}} = - \frac{dV^{\prime}}{dr^{\prime}}$$

To complete the choice of dimensionless variables, we define as unit or time $\tau =\sqrt{m \sigma^{2}/\epsilon}$. This definition leads to a dimensionless form of the equation of motion.

In the boundary condition, we need the ratio $mv^{2}/k_{B}T$. From the dimensionless units of space and time, a dimensionless velocity follows $v = (\sigma/\tau) v^{\prime}$. Therefore we find 
$$ \frac{mv^{2}}{k_{B}T} = \frac{m (\sigma/\tau)^{2} v^{\prime \; 2}}{k_{B}T} = \frac{\epsilon}{k_{B} T} v^{\prime}$$ If we further define the temperature scale in units of $k_{B}T$: $\epsilon^{\prime} = \epsilon/(k_{B}T)$, then we have a complete set of dimensionless equations.

I would choose $L^{\prime} = 10$, $\epsilon^{\prime} = 1$ and about $N = 100$ particles.
{{justifyright{
''Due Friday, November 18, 2016''
}}}

!!! 1. Wave equation in a string of nonuniform density and tension
Let $y(x,t)$ be the vertical displacement of an elastic string from the horizontal, as a function of the lateral coordinate $x$, and time $t$. If the density of the string is a function of position $\rho = \rho(x)$, and the local tension on the string is also a function of position $T=T(x)$, then the linear (small displacements) equation for the string displacement is
$$
\frac{d T(x)}{dx} \frac{\partial y(x,t)}{\partial x} + T(x) \frac{\partial^{2} y(x,t)}{\partial x^{2}} = \rho(x)  \frac{\partial^{2} y(x,t)}{\partial t^{2}},
$$ 
where both $\rho(x) = \rho_{0}e^{\alpha x}$ and $T(x)=T_{0}e^{\alpha x}$ are assumed given.

Write a code that will integrate this equation in $x \in (0,1)$ with $y(0,t)=y(1,t) = 0$. Choose some initial condition that will yield interesting behavior, and plot the resulting evolution.

!!! 2. Wave packet in a quantum well
The time dependent Schroedinger equation for the complex wave function $\psi(x,t)$ is ($\hbar^{2}/2m = 1$),
$$
i \frac{\partial \psi(x,t)}{\partial t} = - \frac{\partial^{2} \psi(x,t)}{\partial x^{2}} + V(x) \psi(x,t)
$$
Consider an electron in  a well which is initially localized at $x=5$ moving with momentum $k_{0} $,
$$
\Large{\psi(x,t=0) = e^{-\frac{1}{2} \left( \frac{x-5}{\sigma_{0}} \right)^{2}} e^{i k_{0} x}}
$$

# Write a code that solves for the motion of the electron in a square  well with$V=0$ in $x\in (0,15)$, with boundary conditions $\psi(0,t) = \psi(15,t) = 0$. Take $\sigma_{0} = 0.5$, $\Delta x = 0.02$, $k_{0} = 17 \pi$, and $\Delta t=\Delta x^{2}/2$. Remember that the equation is complex, so you will have to start by writing its real and imaginary parts first, and then integrate the coupled equations that result.
# Is $\psi(x,t=0)$ an eigenstate of the Hamiltonian?
# Plot the probability distribution that represents the motion of the electron as a function of time.
{{justifyright{
''Due Friday, November 13, 2018''
}}}

!!! 1. Wave equation in a string of nonuniform density and tension
Let $y(x,t)$ be the vertical displacement of an elastic string from the horizontal, as a function of the lateral coordinate $x$, and time $t$. If the density of the string is a function of position $\rho = \rho(x)$, and the local tension on the string is also a function of position $T=T(x)$, then the linear (small displacements) equation for the string displacement is
$$
\frac{d T(x)}{dx} \frac{\partial y(x,t)}{\partial x} + T(x) \frac{\partial^{2} y(x,t)}{\partial x^{2}} = \rho(x)  \frac{\partial^{2} y(x,t)}{\partial t^{2}},
$$ 
where both $\rho(x) = \rho_{0}e^{\alpha x}$ and $T(x)=T_{0}e^{\alpha x}$ are assumed given.

Write a code that will integrate this equation in $x \in (0,1)$ with $y(0,t)=y(1,t) = 0$. Choose some initial condition that will yield interesting behavior, and plot the resulting evolution.

!!! 2. Wave packet in a quantum well
The time dependent Schroedinger equation for the complex wave function $\psi(x,t)$ is ($\hbar^{2}/2m = 1$),
$$
i \frac{\partial \psi(x,t)}{\partial t} = - \frac{\partial^{2} \psi(x,t)}{\partial x^{2}} + V(x) \psi(x,t)
$$
Consider an electron in  a well which is initially localized at $x=5$ moving with momentum $k_{0} $,
$$
\Large{\psi(x,t=0) = e^{-\frac{1}{2} \left( \frac{x-5}{\sigma_{0}} \right)^{2}} e^{i k_{0} x}}
$$

# Write a code that solves for the motion of the electron in a square  well with$V=0$ in $x\in (0,15)$, with boundary conditions $\psi(0,t) = \psi(15,t) = 0$. Take $\sigma_{0} = 0.5$, $\Delta x = 0.02$, $k_{0} = 17 \pi$, and $\Delta t=\Delta x^{2}/2$. Remember that the equation is complex, so you will have to start by writing its real and imaginary parts first, and then integrate the coupled equations that result.
# Is $\psi(x,t=0)$ an eigenstate of the Hamiltonian?
# Plot the probability distribution that represents the motion of the electron as a function of time.
{{justifyright{
''Due Tuesday, December 9, 2014''
}}}
!!Monte Carlo simulation
Develop a Monte Carlo algorithm to simulate the equilibrium distribution of the Ising model in two dimensions. Consider a square lattice with $N \times N$ sites, in which each lattice site $(i,j)$ is occupied by a spin variable adopting the values $S{ij} = \pm 1$.
* The energy of a given configuration $\{ S_{ij} \}$ is given by, $$ E(\{ S_{ij} \}) = -J \sum_{<n.n.>} S_{ij}S_{kl}$$
where the constant $J$ is the interaction constant, and the sum extends only over nearest neighbors in the square lattice, without double counting interection bonds (there are four nearest neighbors in the square lattice: north, south, east and west of the spin in question).
* Use dimensionless variables: $J^{\prime} = J/k_{B}T$.
* Use periodic boundary conditions in both spatial directions. The following modular arithmetic may be useful. If $i$ is a given index (in either direction), then $ip1=i+1$ can be expressed as:$$ip1 = 1 + MOD(i,N)$$.
* Starting from a random initial condition, run two simulations for $J^{\prime} = 0.4, 0.5$. Use the Metropolis algorithm as discussed in class.
* Compute the average energy per spin $\langle E \rangle/N^{2}$ and the total magnetization per spin $\langle M \rangle = (1/N^{2}) \langle \sum_{i,j} S_{i,j} \rangle$.
{{justifyright{
''Due Friday, December 2, 2016''
}}}

!!Molecular Dynamics simulation
Write a Molecular Dynamics simulation code with the following features:
* Consider a fluid system in a ''two dimensional'' square box of size $L \times L$.
* Particles interact via a standard ~Lennard-Jones potential $$V(r) = 4 \epsilon \left( \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r}\right)^{6} \right)$$
where $\epsilon$ and $\sigma$ are constants and assumed known.
* Use the velocity-Verlet algorithm for time integration.
* The system will be held at constant temperature $T$. In order to accomplish that, special boundary conditions will be assumed at the walls of the box (a thermostat at constant temperature). When a particle collides with any of the walls, the tangential component of the velocity of the particle is conserved, and the incident normal velocity $v_{n}$ is updated into an outgoing normal velocity $v_{n}^{\prime}$ given by [ [[van Beijeren|http://arxiv.org/abs/1411.2983]], 2014]: $$v_{n}^{\prime} = \sqrt{ - \frac{2}{\beta m} \ln \left( 1 - {\rm exp} \left( - \frac{\beta m v_{n}^{2}}{2} \right) \right) }$$ Or course, $\beta = 1 /k_{B}T$, and $m$ is the mass of the particles, all assumed identical.
* Make sure you use appropriate dimensionless variables.

Run the code and compute the internal energy, pressure and pair correlation function at some reasonable temperature.

!!!Notes on dimensionless units
Choose $\sigma$ as the unit of length: $r^{\prime}=r/\sigma$ is the dimensionless distance, and $L^{\prime}=L/\sigma$ the dimensionless lateral dimension of the box (e.g. 10). If we introduce a dimensionless potential energy $V^{\prime}=V/\epsilon$ and a dimensionless time $t^{\prime}=t/\tau$, then Newton's Law transforms to
$$ \frac{m \sigma^{2}}{\tau^{2} \epsilon} \frac{d^{2} r^{\prime}}{dt^{\prime \; 2}} = - \frac{dV^{\prime}}{dr^{\prime}}$$

To complete the choice of dimensionless variables, we define as unit or time $\tau =\sqrt{m \sigma^{2}/\epsilon}$. This definition leads to a dimensionless form of the equation of motion.

In the boundary condition, we need the ratio $mv^{2}/k_{B}T$. From the dimensionless units of space and time, a dimensionless velocity follows $v = (\sigma/\tau) v^{\prime}$. Therefore we find 
$$ \frac{mv^{2}}{k_{B}T} = \frac{m (\sigma/\tau)^{2} v^{\prime \; 2}}{k_{B}T} = \frac{\epsilon}{k_{B} T} v^{\prime}$$ If we further define the temperature scale in units of $k_{B}T$: $\epsilon^{\prime} = \epsilon/(k_{B}T)$, then we have a complete set of dimensionless equations.

I would choose $L^{\prime} = 10$, $\epsilon^{\prime} = 1$ and about $N = 100$ particles.
{{justifyright{
''Due Friday, December 9, 2016''
}}}
''Note'': You are not required to submit this homework. If you do by the due date given above, it will be graded and given back to you. The final average homework grade for this class will be computed from your best eight grades. 

!!Monte Carlo simulation
Develop a Monte Carlo algorithm to simulate the equilibrium distribution of the Ising model in two dimensions. Consider a square lattice with $N \times N$ sites, in which each lattice site $(i,j)$ is occupied by a spin variable adopting the values $S{ij} = \pm 1$.
* The energy of a given configuration $\{ S_{ij} \}$ is given by, $$ E(\{ S_{ij} \}) = -J \sum_{<n.n.>} S_{ij}S_{kl}$$
where the constant $J$ is the interaction constant, and the sum extends only over nearest neighbors in the square lattice, without double counting interection bonds (there are four nearest neighbors in the square lattice: north, south, east and west of the spin in question).
* Use dimensionless variables: $J^{\prime} = J/k_{B}T$.
* Use periodic boundary conditions in both spatial directions. The following modular arithmetic may be useful. If $i$ is a given index (in either direction), then the next index, $ip1$ (i plus one), and the prior index $im1$ (i minus 1) can be expressed as:$$ip1 = 1 + MOD(i,N) \quad\quad im1 = 1+MOD(i+N-2,N)$$, where $MOD$ is the MOD function (MOD(a,b) in FORTRAN, a%b in C).
* Starting from a random initial condition, run two simulations for $J^{\prime} = 0.4, 0.5$. Use the Metropolis algorithm as discussed in class.
* Compute the average energy per spin $\langle E \rangle/N^{2}$ and the total magnetization per spin $\langle M \rangle = (1/N^{2}) \langle \sum_{i,j} S_{i,j} \rangle$.
!!!!Homeworks
* [[Homework 1|Homework 1 2018]]. Due on Tuesday, September 18, 2018.
* [[Homework 2|Homework 2 2018]]. Due on Tuesday, September 25, 2018.
* [[Homework 3|Homework 3 2018]]. Due on Tuesday, October 2, 2018.
* [[Homework 4|Homework 4 2018]]. Due on Tuesday, October 9, 2018.
* [[Homework 5|Homework 5 2018]]. Due on Tuesday, October 30,2018.
* [[Homework 6|Homework 6 2018]]. Due on Tuesday, November 6, 2018.
* [[Homework 7|Homework 7 2018]]. Due on Tuesday, November 13, 2018.

!!!!Laboratory sessions
* [[Laboratory Session 1|LabSession 1 2018]]. September 6, 2018.
* [[Laboratory Session 2|LabSession 2 2018]]. September 13, 2018.
* [[Laboratory Session 3|LabSession 3 2018]]. September 20,2018.
* [[Laboratory Session 4|LabSession 4 2018]]. September 27, 2018.
* [[Laboratory Session 5|LabSession 5 2018]]. October 4, 2018.
* [[Laboratory Session 6|LabSession 6 2018]]. October 11, 2018. 
* [[Laboratory Session 7|LabSession 7 2018]]. October 18, 2018.
* [[Laboratory Session 8|LabSession 8 2018]]. October 25, 2018.
* [[Laboratory Session 9|LabSession 9 2018]]. November 1, 2018.
* [[Laboratory Session 10|LabSession 10 2018]]. November 8, 2018.

!!!! Mid term exam
* Take home [[exam |Midterm Exam 2018]] due on Friday, October 19, 2018, 5:00 p.m.
/***
|Name|ImageSizePlugin|
|Source|http://www.TiddlyTools.com/#ImageSizePlugin|
|Version|1.2.3|
|Author|Eric Shulman|
|License|http://www.TiddlyTools.com/#LegalStatements|
|~CoreVersion|2.1|
|Type|plugin|
|Description|adds support for resizing images|
This plugin adds optional syntax to scale an image to a specified width and height and/or interactively resize the image with the mouse.
!!!!!Usage
<<<
The extended image syntax is:
{{{
[img(w+,h+)[...][...]]
}}}
where ''(w,h)'' indicates the desired width and height (in CSS units, e.g., px, em, cm, in, or %). Use ''auto'' (or a blank value) for either dimension to scale that dimension proportionally (i.e., maintain the aspect ratio). You can also calculate a CSS value 'on-the-fly' by using a //javascript expression// enclosed between """{{""" and """}}""". Appending a plus sign (+) to a dimension enables interactive resizing in that dimension (by dragging the mouse inside the image). Use ~SHIFT-click to show the full-sized (un-scaled) image. Use ~CTRL-click to restore the starting size (either scaled or full-sized).
<<<
!!!!!Examples
<<<
{{{
[img(100px+,75px+)[images/meow2.jpg]]
}}}
[img(100px+,75px+)[images/meow2.jpg]]
{{{
[<img(34%+,+)[images/meow.gif]]
[<img(21% ,+)[images/meow.gif]]
[<img(13%+, )[images/meow.gif]]
[<img( 8%+, )[images/meow.gif]]
[<img( 5% , )[images/meow.gif]]
[<img( 3% , )[images/meow.gif]]
[<img( 2% , )[images/meow.gif]]
[img(  1%+,+)[images/meow.gif]]
}}}
[<img(34%+,+)[images/meow.gif]]
[<img(21% ,+)[images/meow.gif]]
[<img(13%+, )[images/meow.gif]]
[<img( 8%+, )[images/meow.gif]]
[<img( 5% , )[images/meow.gif]]
[<img( 3% , )[images/meow.gif]]
[<img( 2% , )[images/meow.gif]]
[img(  1%+,+)[images/meow.gif]]
{{tagClear{
}}}
<<<
!!!!!Revisions
<<<
2011.09.03 [1.2.3] bypass addStretchHandlers() if no '+' suffix is used (i.e., not resizable)
2010.07.24 [1.2.2] moved tip/dragtip text to config.formatterHelpers.imageSize object to enable customization
2009.02.24 [1.2.1] cleanup width/height regexp, use '+' suffix for resizing
2009.02.22 [1.2.0] added stretchable images
2008.01.19 [1.1.0] added evaluated width/height values
2008.01.18 [1.0.1] regexp for "(width,height)" now passes all CSS values to browser for validation
2008.01.17 [1.0.0] initial release
<<<
!!!!!Code
***/
//{{{
version.extensions.ImageSizePlugin= {major: 1, minor: 2, revision: 3, date: new Date(2011,9,3)};
//}}}
//{{{
var f=config.formatters[config.formatters.findByField("name","image")];
f.match="\\[[<>]?[Ii][Mm][Gg](?:\\([^,]*,[^\\)]*\\))?\\[";
f.lookaheadRegExp=/\[([<]?)(>?)[Ii][Mm][Gg](?:\(([^,]*),([^\)]*)\))?\[(?:([^\|\]]+)\|)?([^\[\]\|]+)\](?:\[([^\]]*)\])?\]/mg;
f.handler=function(w) {
	this.lookaheadRegExp.lastIndex = w.matchStart;
	var lookaheadMatch = this.lookaheadRegExp.exec(w.source)
	if(lookaheadMatch && lookaheadMatch.index == w.matchStart) {
		var floatLeft=lookaheadMatch[1];
		var floatRight=lookaheadMatch[2];
		var width=lookaheadMatch[3];
		var height=lookaheadMatch[4];
		var tooltip=lookaheadMatch[5];
		var src=lookaheadMatch[6];
		var link=lookaheadMatch[7];

		// Simple bracketted link
		var e = w.output;
		if(link) { // LINKED IMAGE
			if (config.formatterHelpers.isExternalLink(link)) {
				if (config.macros.attach && config.macros.attach.isAttachment(link)) {
					// see [[AttachFilePluginFormatters]]
					e = createExternalLink(w.output,link);
					e.href=config.macros.attach.getAttachment(link);
					e.title = config.macros.attach.linkTooltip + link;
				} else
					e = createExternalLink(w.output,link);
			} else 
				e = createTiddlyLink(w.output,link,false,null,w.isStatic);
			addClass(e,"imageLink");
		}

		var img = createTiddlyElement(e,"img");
		if(floatLeft) img.align="left"; else if(floatRight) img.align="right";
		if(width||height) {
			var x=width.trim(); var y=height.trim();
			var stretchW=(x.substr(x.length-1,1)=='+'); if (stretchW) x=x.substr(0,x.length-1);
			var stretchH=(y.substr(y.length-1,1)=='+'); if (stretchH) y=y.substr(0,y.length-1);
			if (x.substr(0,2)=="{{")
				{ try{x=eval(x.substr(2,x.length-4))} catch(e){displayMessage(e.description||e.toString())} }
			if (y.substr(0,2)=="{{")
				{ try{y=eval(y.substr(2,y.length-4))} catch(e){displayMessage(e.description||e.toString())} }
			img.style.width=x.trim(); img.style.height=y.trim();
			if (stretchW||stretchH) config.formatterHelpers.addStretchHandlers(img,stretchW,stretchH);
		}
		if(tooltip) img.title = tooltip;

		// GET IMAGE SOURCE
		if (config.macros.attach && config.macros.attach.isAttachment(src))
			src=config.macros.attach.getAttachment(src); // see [[AttachFilePluginFormatters]]
		else if (config.formatterHelpers.resolvePath) { // see [[ImagePathPlugin]]
			if (config.browser.isIE || config.browser.isSafari) {
				img.onerror=(function(){
					this.src=config.formatterHelpers.resolvePath(this.src,false);
					return false;
				});
			} else
				src=config.formatterHelpers.resolvePath(src,true);
		}
		img.src=src;
		w.nextMatch = this.lookaheadRegExp.lastIndex;
	}
}

config.formatterHelpers.imageSize={
	tip: 'SHIFT-CLICK=show full size, CTRL-CLICK=restore initial size',
	dragtip: 'DRAG=stretch/shrink, '
}

config.formatterHelpers.addStretchHandlers=function(e,stretchW,stretchH) {
	e.title=((stretchW||stretchH)?this.imageSize.dragtip:'')+this.imageSize.tip;
	e.statusMsg='width=%0, height=%1';
	e.style.cursor='move';
	e.originalW=e.style.width;
	e.originalH=e.style.height;
	e.minW=Math.max(e.offsetWidth/20,10);
	e.minH=Math.max(e.offsetHeight/20,10);
	e.stretchW=stretchW;
	e.stretchH=stretchH;
	e.onmousedown=function(ev) { var ev=ev||window.event;
		this.sizing=true;
		this.startX=!config.browser.isIE?ev.pageX:(ev.clientX+findScrollX());
		this.startY=!config.browser.isIE?ev.pageY:(ev.clientY+findScrollY());
		this.startW=this.offsetWidth;
		this.startH=this.offsetHeight;
		return false;
	};
	e.onmousemove=function(ev) { var ev=ev||window.event;
		if (this.sizing) {
			var s=this.style;
			var currX=!config.browser.isIE?ev.pageX:(ev.clientX+findScrollX());
			var currY=!config.browser.isIE?ev.pageY:(ev.clientY+findScrollY());
			var newW=(currX-this.offsetLeft)/(this.startX-this.offsetLeft)*this.startW;
			var newH=(currY-this.offsetTop )/(this.startY-this.offsetTop )*this.startH;
			if (this.stretchW) s.width =Math.floor(Math.max(newW,this.minW))+'px';
			if (this.stretchH) s.height=Math.floor(Math.max(newH,this.minH))+'px';
			clearMessage(); displayMessage(this.statusMsg.format([s.width,s.height]));
		}
		return false;
	};
	e.onmouseup=function(ev) { var ev=ev||window.event;
		if (ev.shiftKey) { this.style.width=this.style.height=''; }
		if (ev.ctrlKey)  { this.style.width=this.originalW; this.style.height=this.originalH; }
		this.sizing=false;
		clearMessage();
		return false;
	};
	e.onmouseout=function(ev) { var ev=ev||window.event;
		this.sizing=false;
		clearMessage();
		return false;
	};
}
//}}}
[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/Kmeans1d.jpg]]

[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/Kmeans2d.jpg]]
#  Try your login credentials, or create them if you have not done so.
# Write a simple code, compile it an run it that will estimate the largest simple precision (4 byte) and double precision (8 byte) floating point and integer numbers that the computer you are using can represent. Hint: Use a loop to compute 1*2*2*2 .... until the returned value is out of range.
# The Bessel function of order 1 is given by $$ J_{1}(x) = \frac{\sin x}{x^{2}} - \frac{\cos x}{x}$$ Given the Taylor series for sin and cos, $$ \sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - ... \quad \cos (x) = 1 - \frac{x^{2}}{2!} + \frac{x^{4}}{4!} - ...$$Write a code that computes this function in $x \in [0,10]$ and plot the result.
# Try your x500 login credentials. Ideally you should have tried them before you come to the lab.
# Log in github.umn.edu with your x500 credentials. Inform the instructor by e-mail when you have done so. 
# Write a simple code, compile it an run it that will estimate the largest simple precision (4 byte) and double precision (8 byte) floating point and integer numbers that the computer you are using can represent. Hint: Use a loop to compute 1*2*2*2 .... until the returned value is out of range.
# The Bessel function of order 1 is given by $$ J_{1}(x) = \frac{\sin x}{x^{2}} - \frac{\cos x}{x}$$ Given the Taylor series for sin and cos, $$ \sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - ... \quad \cos (x) = 1 - \frac{x^{2}}{2!} + \frac{x^{4}}{4!} - ...$$Write a code that uses the Taylor series to compute the Bessel function in $x \in [0,10]$ to a reasonable approximation and plot the result.
{{justifyright{
''Thursday, September 6, 2018''
}}}

# Try your x500 login credentials. Ideally you should have tried them before you come to the lab.
# Log in github.umn.edu with your x500 credentials. Inform the instructor by e-mail when you have done so. 
# Write a simple code, compile it an run it that will estimate the largest simple precision (4 byte) and double precision (8 byte) floating point and integer numbers that the computer you are using can represent. Hint: Use a loop to compute 1*2*2*2 .... until the returned value is out of range.
# The Bessel function of order 1 is given by $$ J_{1}(x) = \frac{\sin x}{x^{2}} - \frac{\cos x}{x}$$ Given the Taylor series for sin and cos, $$ \sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - ... \quad \cos (x) = 1 - \frac{x^{2}}{2!} + \frac{x^{4}}{4!} - ...$$Write a code that uses the Taylor series to compute the Bessel function in $x \in [0,10]$ to a reasonable approximation and plot the result.

Compilation instructions:
* ''FORTRAN'' (type man gfortran for information and options. Information [[about Fortran file extensions|http://fortranwiki.org/fortran/show/File+extensions]]).
** gfortran filename.f90
* ''C/C++'' (type man g++ or man gcc for information and options)
** g++ filename.c or gcc filename.c -lm
!!Random numbers and their distributions
So called pseudo-random numbers can be generated on a computer. There are multiple routines that can accomplih this to varying degrees of sophistication. We will use a simple intrinsic function on Linux computers. The function call (both FORTRAN and C)
{{{
x=RAND()
}}}
assigns to x a random number uniformly distributed between $[0,1)$ (any number in that interval with equal probability).

Because the numbers are not truly random (rather they are a predictable sequence), a sequence has to be initialized with an integer number called 'seed'. For example
{{{
INTEGER SEED
SEED = 123456
CALL SRAND(SEED)
...
x=RAND()
}}}
Every time the same seed is used, the same sequence of random numbers results. A sequence of random numbers is obtained simply by repeating the function call RAND().

!!!!1. Uniform random number statistics
Just to get started, generate 100,000 random numbers and compute the average and standard deviation of the sample. You should be able to predict the values analytically and compare with the results.

!!!!2. Histogram
We will create a histogram of the sequence of random numbers and compare it to the expected probability (uniform). In order to create a histogram, you can use the following procedure
{{{
INTEGER :: hist(100)
REAL(8) :: rhist(100)

DO i=1,100
hist(i)=0
END DO

dx=(1.0-0.0)/100
CALL SRAND(3452345)

DO i=1,100000
index=RAND()/dx
hist(index)=hist(index)+1
END DO

DO i=1,100
rhist(i)=hist(i)/(dx*100000)
WRITE(27,*) i*dx,rhist(i)
END DO
}}}
Plot the result and compare the histogram to the expected probability distribution of the random numbers.

!!!!3. Gaussianly distributed random numbers
Gaussianly distributed random numbers are generated in pairs (we will see this in class). Given two random numbers $r_{1}$ and $r_{2}$ uniformly distributed in [0,1), the following two numbers
$$x_{1} = \sqrt{-2 \ln r_{1}} \; \cos (2 \pi r_{2}) \quad\quad x_{2} = \sqrt{-2 \ln r_{1}} \; \sin (2 \pi r_{2})$$
are two independent random numbers, Gaussianly distributed ov zero average and variance equal to 1.

Repeat the process of (2) to generate the histogram of the resulting sequence.

!!!!4. Correlations
The sequences above involve independent random numbers and therefore there are no correlations. We will examine this first. Assume that you have generated the sequence of random numbers $\{ x_{i} \}$. The correlation between, say, $x_{i}$ and $x_{i+1}$ is defined
$$C_{i\;i+1} = \langle x_{i} x_{i+1} \rangle - \langle x_{i} \rangle \langle x_{i+1} \rangle$$
where average over $N$ points are defined as
$$\langle ... \rangle \approx \frac{1}{N} \sum_{i=1}^{N} (...)_{i}$$

Calculate $C_{i \; i+1}$ for the Gaussian sequence above.

Correlation functions are a generalization in that one defines
$$C(j) = \langle x_{i} x_{i+j} \rangle - \langle x_{i} \rangle \langle x_{i+j} \rangle$$ 
Note that the correlation only depends on $j$, the "distance" between the two points. This is normally the case in Physics, but it is not always the case.

!!!!5. More correlations
Consider the following sequence
$$x_{i+1} = x_{i} - x_{i-1} + r_{i}$$
where $r_{i}$ is a random number uniformly distributed in [0,1). Specify some initial condition $x_{0}$ and $x_{1}$. Run a long sequence and calculate the correlation function $C(j)$ for $j$ between 1 and 10.
* [[Presentation|http://homepages.spa.umn.edu/~vinals/tspot_files/Lab10.pdf]] for this session posted.
* We will next cover some [[Programming Tips]].
* The final code has been added to github.umn.edu.
{{justifyright{
''Thursday, November 8, 2018''
}}}
!!Random numbers and their distributions
So called pseudo-random numbers can be generated on a computer. There are multiple routines that can accomplish this to varying degrees of sophistication. We will use a simple intrinsic function on Linux computers. A function call (both in FORTRAN and C++) assigns a random number uniformly distributed between $[0,1)$ (any number in that interval with equal probability). Details on the syntax can be found [[here|http://homepages.spa.umn.edu/~vinals/classes/phys4041/phys4041_notes.pdf]].

Go over the sequence of steps indicated below as time allows during the laboratory session.

!!!!1. Uniform random number statistics
Just to get started, generate 100,000 random numbers and compute the average and standard deviation of the sample. You should be able to predict the values analytically and compare with the results.

!!!!2. Histogram
We will create a histogram of the sequence of random numbers and compare it to the expected probability (uniform). In order to create a histogram, you can use the following procedure (e.g., in FORTRAN. Similarly in C++)
{{{
INTEGER :: hist(100)
REAL(8) :: rhist(100)

DO i=1,100
hist(i)=0
END DO

dx=(1.0-0.0)/100
CALL SRAND(3452345)

DO i=1,100000
index=RAND()/dx + 1
hist(index)=hist(index)+1
END DO

DO i=1,100
rhist(i)=hist(i)/(dx*100000)
WRITE(27,*) i*dx,rhist(i)
END DO
}}}
Plot the result and compare the histogram to the expected probability distribution of the random numbers.

!!!!3. Gaussianly distributed random numbers
Gaussianly distributed random numbers are generated in pairs (we will see this in class). Given two random numbers $r_{1}$ and $r_{2}$ uniformly distributed in [0,1), the following two numbers
$$x_{1} = \sqrt{-2 \ln r_{1}} \; \cos (2 \pi r_{2}) \quad\quad x_{2} = \sqrt{-2 \ln r_{1}} \; \sin (2 \pi r_{2})$$
are two independent random numbers, Gaussianly distributed of zero average and variance equal to 1.

Repeat the process of (2) to generate the histogram of the resulting sequence.

!!!!4. Correlations
The sequences above involve independent random numbers and therefore there are no correlations. We will examine this first. Assume that you have generated the sequence of random numbers $\{ x_{i} \}$. The correlation between, say, $x_{i}$ and $x_{i+1}$ is defined
$$C_{i\;i+1} = \langle x_{i} x_{i+1} \rangle - \langle x_{i} \rangle \langle x_{i+1} \rangle$$
where average over $N$ points are defined as
$$\langle ... \rangle \approx \frac{1}{N} \sum_{i=1}^{N} (...)_{i}$$

Calculate $C_{i \; i+1}$ for the Gaussian sequence above.

Correlation functions are a generalization in that one defines
$$C(j) = \langle x_{i} x_{i+j} \rangle - \langle x_{i} \rangle \langle x_{i+j} \rangle$$ 
Note that the correlation only depends on $j$, the "distance" between the two points. This is normally the case in Physics, but it is not always the case.

!!!!5. More correlations
Consider the following sequence
$$x_{i+1} = x_{i} - x_{i-1} + r_{i}$$
where $r_{i}$ is a random number uniformly distributed in [0,1). Specify some initial condition $x_{0}$ and $x_{1}$. Run a long sequence and calculate the correlation function $C(j)$ for $j$ between 1 and 10.
!!!Correlation Functions
!!!!!Independent and Gaussianly distributed random numbers
Gaussianly distributed random numbers are generated in pairs (as seen in class). Given two random numbers $r_{1}$ and $r_{2}$ uniformly distributed in [0,1), the following two numbers
$$x_{1} = \sqrt{-2 \ln r_{1}} \; \cos (2 \pi r_{2}) \quad\quad x_{2} = \sqrt{-2 \ln r_{1}} \; \sin (2 \pi r_{2})$$
are two independent random numbers, Gaussianly distributed ov zero average and variance equal to 1. We wish to check both the Gaussianity and the independence.

In order to check for Gaussianity, compute $\langle x \rangle$, $\langle x^{2} \rangle$,$\langle x^{3} \rangle$ and $\langle x^{4} \rangle$ where the averages are defined over the sequence:
$$\langle ... \rangle \approx \frac{1}{N} \sum_{i=1}^{N} (...)_{i}$$
Show that odd moments are zero, that $\langle x^{2} \rangle = \sigma^{2} = 1$ (the variance of the sequence), and that $\langle x^{4} \rangle = 3 \langle x^{2} \rangle$ as appropriate for a Gaussian distribution.

In order to check for independence in a sequence of random numbers $\{ x_{i} \}$, we compute the correlation between, say, $x_{i}$ and $x_{i+1}$ defined by
$$C(1) = \langle x_{i} x_{i+1} \rangle - \langle x_{i} \rangle \langle x_{i+1} \rangle$$
Or more generally, we can define a correlation function
$$C(j) = \langle x_{i} x_{i+j} \rangle - \langle x_{i} \rangle \langle x_{i+j} \rangle$$

Calculate $C(j)$, $j=1, ..., 10$ for the sequence of Gaussian numbers above.

!!!!!Correlations in an oscillatory process
In general, stochastic processes generate states which are correlated. As an example, consider the sequence,
$$x_{i+1} = x_{i} - x_{i-1} + r_{i}$$
where $r_{i}$ is a sequence of independent random numbers uniformly distributed in [0,1). Specify some initial condition $x_{0}$ and $x_{1}$. This iteration produces a sequence of random numbers $\{ x_{i} \}$.

Run a long sequence and plot it. Next, calculate the correlation function $C(j)$ for $j$ between 1 and 10.
!!Random numbers and their distributions
So called pseudo-random numbers can be generated on a computer. There are multiple routines that can accomplih this to varying degrees of sophistication. We will use a simple intrinsic function on Linux computers. The function call (both FORTRAN and C)
{{{
x=RAND()
}}}
assigns to x a random number uniformly distributed between $[0,1)$ (any number in that interval with equal probability).

Because the numbers are not truly random (rather they are a predictable sequence), a sequence has to be initialized with an integer number called 'seed'. For example
{{{
INTEGER SEED
SEED = 123456
CALL SRAND(SEED)
...
x=RAND()
}}}
Every time the same seed is used, the same sequence of random numbers results. A sequence of random numbers is obtained simply by repeating the function call RAND().

Go over the sequence steps below as time allows during the laboratory session.

!!!!1. Uniform random number statistics
Just to get started, generate 100,000 random numbers and compute the average and standard deviation of the sample. You should be able to predict the values analytically and compare with the results.

!!!!2. Histogram
We will create a histogram of the sequence of random numbers and compare it to the expected probability (uniform). In order to create a histogram, you can use the following procedure
{{{
INTEGER :: hist(100)
REAL(8) :: rhist(100)

DO i=1,100
hist(i)=0
END DO

dx=(1.0-0.0)/100
CALL SRAND(3452345)

DO i=1,100000
index=RAND()/dx
hist(index)=hist(index)+1
END DO

DO i=1,100
rhist(i)=hist(i)/(dx*100000)
WRITE(27,*) i*dx,rhist(i)
END DO
}}}
Plot the result and compare the histogram to the expected probability distribution of the random numbers.

!!!!3. Gaussianly distributed random numbers
Gaussianly distributed random numbers are generated in pairs (we will see this in class). Given two random numbers $r_{1}$ and $r_{2}$ uniformly distributed in [0,1), the following two numbers
$$x_{1} = \sqrt{-2 \ln r_{1}} \; \cos (2 \pi r_{2}) \quad\quad x_{2} = \sqrt{-2 \ln r_{1}} \; \sin (2 \pi r_{2})$$
are two independent random numbers, Gaussianly distributed of zero average and variance equal to 1.

Repeat the process of (2) to generate the histogram of the resulting sequence.

!!!!4. Correlations
The sequences above involve independent random numbers and therefore there are no correlations. We will examine this first. Assume that you have generated the sequence of random numbers $\{ x_{i} \}$. The correlation between, say, $x_{i}$ and $x_{i+1}$ is defined
$$C_{i\;i+1} = \langle x_{i} x_{i+1} \rangle - \langle x_{i} \rangle \langle x_{i+1} \rangle$$
where average over $N$ points are defined as
$$\langle ... \rangle \approx \frac{1}{N} \sum_{i=1}^{N} (...)_{i}$$

Calculate $C_{i \; i+1}$ for the Gaussian sequence above.

Correlation functions are a generalization in that one defines
$$C(j) = \langle x_{i} x_{i+j} \rangle - \langle x_{i} \rangle \langle x_{i+j} \rangle$$ 
Note that the correlation only depends on $j$, the "distance" between the two points. This is normally the case in Physics, but it is not always the case.

!!!!5. More correlations
Consider the following sequence
$$x_{i+1} = x_{i} - x_{i-1} + r_{i}$$
where $r_{i}$ is a random number uniformly distributed in [0,1). Specify some initial condition $x_{0}$ and $x_{1}$. Run a long sequence and calculate the correlation function $C(j)$ for $j$ between 1 and 10.
Continue debugging your Molecular Dynamics code.

Once it runs, experiment with:
* A small number of particles $N = 10$ in a small box $L=10$. Plot trajectories and examine the effect of boundary conditions. Hint: an easy way to plot trajectories is to write to a file with $(x,y)$ for all particles and a sequence of times and simply use gnuplot.
* Compare the trajectories between the boundary conditions given and a different condition in which particles are just reflected off boundaries.
* Examine the effect of temperature on the trajectories, as well as the values of the energy and pressure.
!!! Correlation function of the two dimensional ~Lennard-Jones system
Use the Molecular Dynamics code that you have written for the homework to compute the pairwise correlation function of the ~Lennard-Jones system for a range of temperatures. The correlation function is defined as a histogram of inter particle separations $$ g(\vec{r}) = \frac{V}{N} \langle \sum_{i,j} \delta(\vec{r} - \vec{r}_{ij}) \rangle 
$$
where $\vec{r}_{ij}$ is the vector joining the positions of particles $i$ and $j$, and $j$ is different from $i$. The average is understood over time during the simulation. For simplicity, just focus on separations $r = \| \vec{r} \|$ and only compute the correlation function of separations $g(r)$.

Calculate and plot the correlation function $g(r)$ in the two following states: Fix the dimensionless temperature to $T = 0.45$ and examine two dimensionless densities $\rho = N/V = 0.855$ and $\rho = 0.786$.
!!Gillespie algorithm for a simple production-degradation reaction system
Consider the following two reactions giving production and degradation of X, $$\emptyset \;\; \overset{A}\rightarrow \;\; X \quad\quad X \;\; \overset{B}\rightarrow \;\; \emptyset$$
where production occurs at a constant rate $A$ and degradation at a constant rate $B$.

If $n$ is the number of molecules of $X$, the Master Equation that corresponds to these two reactions is,$$ \frac{\partial p_{n}(t)}{\partial t} = A(E^{-1}-1) p_{n}(t) + B (E-1)(np_{n}(t))$$ 

Identifying the propensitties of the two reactions as $a_{1} = A$ and $a_{2} = Bn$ write a code that implements the Gillespie algorithm to evolve some arbitrary initial condition to the steady state distribution. 

In this case, the steady state distribution is known to be $$ p_{n} = \sqrt{ \frac{B}{2 \pi A}} e^{-\frac{B}{2A} \left( n - \frac{A}{B} \right)^{2}}$$
!Version control and respositories

We will use git as the version control, and the ~UofM github respository. Help on git can be found at https://www.atlassian.com/git

!!Create the respository where class code examples will be stored
* Point your browser to github.umn.edu and login using your X500 credentials.
* Clone the respository on your local machine (you can do this on as many machines as you wish. They will all syncronize). Go the directory in which you would like to have the respository created and type:
{{{
git clone https://github.umn.edu/jvinals/PHYS4041-2016.git
}}}
The subdirectory ~PHYS4041-2016 will be created which will be your local copy of the repository.

Although we will use Linux to conduct class activities, there are also desktop applications available for Mac and windows. You can also synchronize those operating systems if desired. The desktops can be found at  https://desktop.github.com/

You can also create your own repository for your private work. To do so, 
* log in github.umn.edu
* Click on green button "New repository".
* Give it a name (e.g., myrepo), and (typically) choose private.
* Select "Initialize with a README".
* On your local computer, go to the directory where you want the repository to reside and type
{{{
git clone https://github.umn.edu/[your x500 username]/myrepo.git
}}}
A new subdirectory myrepo will be created that can be syncronized with github.umn.edu.
!!Version control
!!!!!0. Synchronize
Synchronization between your local computer and the repository is simple. Once inside the directory that needs to be synchronized, you would type:
{{{
git pull
}}}
to synchronize the local files to those in the repository ("pull" from the repository). 
!!!!!1. Add files that will be versioned to the repository
Use an editor to create the file test.f90. In order to turn on version control, you need to "add" the file to git. Type
{{{
git add test.f90
}}}
git maintains its own list of files subject to revision control. Other files can be created in the same directory that git will simply ignore.

In order to check the status of your files and versions, you can type
{{{
git status
}}}
In the case at hand, it will list a "new file", test.f90
!!!!!2. Committing changes
You can edit and manipulate your file at will. Once you are satisfied with the changes, you can "commit" the latest version to the local repository. git allows you to compare and restore all committed versions in the past. In order to commit, you would type
{{{
git commit
}}}
An editor window will open where you would type any informational message you want kept with the version committed. For example, "do loop fixed and print option added". Any message that will help you later understand what this version is. To finish, hit the Esc key followed by :wq (these are vi commands). Every file in the "add" list that has been changed has now been committed. You can type git status to see where you are.
!!!!!3. Branches
!!!!1. Use of make to maintain scientific code.
Write a makefile that will be used to compile the code you wrote for the homework to check the accuracy of the machine you are using.

Scientific code is usually comprised of a relatively large number of modules and links to external libraries. The command "make" in Linux simplifies code maintenance a great deal. I will list a couple of examples, here. A few more can be found [[here|http://linux.101hacks.com/unix/make/]].

The command make simply executes a file called (by default, this can be changed) makefile. The command simply interprets the statements contained in makefile. The simplest example concerns compiling the code that you produced to estimate the accuracy of the computer you are using. The task of compiling the code would be undertaken by make. Assuming your code was named accuracy.f, create a file makefile with:

{{{
accuracy: accuracy.f
[Tab] gfortran -o accuracy accuracy.f
}}}

Typing make at the command prompt, will compile the code accuracy.f ''only'' when modified.

The general format is:

{{{
Target : Dependencies
[Tab] command to execute if files in dependencies modified
}}}

Target is the final product of the make command. Target depends on other files (e.g. its own source code, or other external codes). If any of the dependencies is satisfied, then the command line following it (after the tab) will execute.

A more complex example in which two codes are part of the same task:
{{{
combined-file: file1.o file2.o
[Tab]   gcc file1.o file2.o -o combined-file

file1.o : file1.c
[Tab]   gcc -c file1.c -o file1.o

file2.o : file2.c
[Tab]    gcc -c file2.c -o file2.o
}}}

make can parse variables, and do more complicated tasks. Here is a more complex example:
{{{
.SUFFIXES = .o .f90 .f
OPTION = OPTIM
COMP_FORTRAN = gfortran
COMP_C = gcc

OBJECTS = file1.o file2.o 

final_code: final_code.f $(OBJECTS)
[Tab]  $(COMP_FORTRAN) $(OPTION) -o final_code final_code.f $(OBJECTS)

.f.o:;
[Tab] $(COMP_FORTRAN) -c $.f
}}}

!!!!2. Git and version control
Modern day coding and code maintenance makes heavy use of version control software and repositories. For this class, we will use git as version control, and github.umn.edu as the repository. For this session we will only use it to download files.
# Open a browser, navigate to github.umn.edu and log on with your x500 credentials.
# @@color(maroon):Contact the instructor to receive access to the repository (after you have logged on at least once)@@.
# Search for the respository jvinals/phys4041.
# Click on the file to download.
# On the top right, click on Raw.
# Save the file from your browser to your computer.

You may also setup a local directory to copy an existing git repository.
# Open a terminal window on your machine. 
# In your home directory (or create/navigate to one less cluttered), enter: {{{ git init }}}
# Enter: {{{ git clone https://github.umn.edu/jvinals/phys4041.git }}}
# You will need to enter your x500 credentials.
# Once downloaded, you will have a local copy of the repository in a directory named phys4041
# To download future updates to the repository, navigate to this directory and enter: {{{ git pull }}}

More info on git can be found in this tutorial: https://www.atlassian.com/git

!!!!3. Linking to external code or libraries
This is a simple test of using external code to accomplish the same task as in the homework. Either:
* From github.umn.edu/jvinals/phys4041 download the subroutines/functions dlamch.f and lsame.f (FORTRAN), dlamch.c, blaswrap.h (C). 
* Write the main code that uses this subroutines/functions to obtain information about the accuracy of the computer. The source code contains help about the various parameters that the functions accept.
Or:
* ''Only on Physics systems'': You can link to the LAPACK package directly that contains those subroutines without any need to dowload. The entire package LAPACK (numerical linear algebra) already compiled is in the archive  /usr/lib64/liblapack.so. If your main code were called dlamch_drive.f90, you would compile by typing:
{{{
gfortran dlamch_drive.f90 -o dlamch_drive /usr/lib64/liblapack.so
}}}
which is often abbreviated into
{{{
gfortran dlamch_drive.f90 -o dlamch_drive -L /usr/lib64/ -llapack
}}}

Create a makefile to automate the code build.
!!!!1. Use of make to maintain scientific code
Write a makefile that will be used to compile the code you wrote during the first Lab session to check the largest and smallest floating point numbers of the machine you are using.

Scientific code is usually comprised of a relatively large number of modules and links to external libraries. The command "make" in Linux simplifies code maintenance a great deal. I will list a couple of examples, here. A few more can be found [[here|http://linux.101hacks.com/unix/make/]].

The command make simply executes a file called (by default, this can be changed) makefile. The command simply interprets the statements contained in makefile. The simplest example concerns compiling the code mentioned above. The task of compiling the code would be undertaken by make. Assuming your code was named accuracy.f, create a file makefile with (FORTRAN)
{{{
accuracy: accuracy.f
[Tab] gfortran -o accuracy accuracy.f
}}}
or C:
{{{
accuracy: accuracy.c
[Tab] g++ -o accuracy accuracy.c
}}}

Typing make at the command prompt, will compile the code accuracy.f or accuracy.c ''only'' when modified, and create the executable named accuracy.

The general format of a  makefile entry is:

{{{
Target: Dependencies
[Tab] command to execute if files in dependencies modified
[Tab] another command to be executed
}}}

Target is the final product of the make command. Target depends on other files (e.g. its own source code, or other external codes). If any of the dependencies is satisfied, then the command line following it (after the tab) will execute.

A more complex example in which two codes are part of the same task:
{{{
combined-file: file1.o file2.o
[Tab]   g++ file1.o file2.o -o combined-file

file1.o: file1.c
[Tab]   g++ -c file1.c -o file1.o

file2.o: file2.c
[Tab]   g++ -c file2.c -o file2.o
}}}

make can parse variables, and do more complicated tasks. Here is a more complex example:
{{{
.SUFFIXES = .o .f90 .f
OPTION = OPTIM
COMP_FORTRAN = gfortran
COMP_C = gcc

OBJECTS = file1.o file2.o 

final_code: final_code.f $(OBJECTS)
[Tab]  $(COMP_FORTRAN) $(OPTION) -o final_code final_code.f $(OBJECTS)

.f.o:;
[Tab] $(COMP_FORTRAN) -c $.f
}}}

!!!!2. Linking to external code or libraries
In this case, you will learn how to use pieces of code (libraries) already written by others. As an example, 
we will use already existing external code to find the accuracy of the machine that you are using.

Either:
* From github.umn.edu/jvinals/phys4041-2018 download the subroutines/functions dlamch.f and lsame.f (FORTRAN), dlamch.c, blaswrap.h (C). {{{git pull}}} should do it if you are set up correctly with github.
* Write the main code that uses this subroutines/functions to obtain information about the accuracy of the computer. The source code of the routines contains help about the various parameters that the functions accept.
Or (the best option):
* Alternatively, you can link to the existing LAPACK package directly that contains those subroutines without any need to download. The entire LAPACK package (numerical linear algebra) already compiled is in the archive  (CSELABS) /usr/lib/liblapack.so for FORTRAN and /usr/lib/liblapacke.so (C), and (Physics) /usr/lib64/liblapack.so (FORTRAN) and /usr/lib64/liblapacke.so (C). For example, if your main code were called dlamch_drive.f90, you would compile it on a Physics department system by typing:
{{{
gfortran dlamch_drive.f90 -o dlamch_drive /usr/lib64/liblapack.so
}}}
which is often abbreviated into
{{{
gfortran dlamch_drive.f90 -o dlamch_drive -L /usr/lib64/ -llapack
}}}
In C, you would type something similar. In the CSELABS systems, you would type
{{{
gcc dlamch_drive.c  -L /usr/lib -llapack
}}}
Create a makefile to automate the code build.
!!!1. Use of make to maintain scientific code
Goal: To write a "makefile" that will be used to compile the code you wrote during the first Lab session to check the largest and smallest floating point numbers of the machine you are using. Scientific code is usually comprised of a relatively large number of modules and links to external libraries. The command "make" in Linux simplifies code maintenance a great deal. I will list a couple of examples here. A few more can be found [[here|http://linux.101hacks.com/unix/make/]].

The command make simply executes a file called (by default, this can be changed) makefile. The command simply interprets the statements contained in makefile. The simplest example concerns compiling the code mentioned above. Assuming your code was named accuracy.f90 or accuracy.c, use a text editor to create a file makefile with the following two lines (FORTRAN)
{{{
accuracy: accuracy.f90
[Tab] gfortran -o accuracy accuracy.f90
}}}
or C:
{{{
accuracy: accuracy.c
[Tab] g++ -o accuracy accuracy.c
}}}

Typing make at the command prompt, will compile the code accuracy.f or accuracy.c ''only'' when modified, and create the executable named accuracy.

The general format of a  makefile entry is:

{{{
Target: Dependencies
[Tab] command to execute if files in dependencies modified
[Tab] another command to be executed
}}}

Target is the final product of the make command. Target depends on other files (e.g. its own source code, or other external codes). If any of the dependencies is satisfied, then the command line following it (after the tab) will execute.

A more complex example in which two codes are part of the same task:
{{{
combined-file: file1.o file2.o
[Tab]   g++ file1.o file2.o -o combined-file

file1.o: file1.c
[Tab]   g++ -c file1.c -o file1.o

file2.o: file2.c
[Tab]   g++ -c file2.c -o file2.o
}}}

make can parse variables, and do more complicated tasks. Here is a more complex example:
{{{
.SUFFIXES = .o .f90 .f
OPTION = OPTIM
COMP_FORTRAN = gfortran
COMP_C = gcc

OBJECTS = file1.o file2.o 

final_code: final_code.f90 $(OBJECTS)
[Tab]  $(COMP_FORTRAN) $(OPTION) -o final_code final_code.f90 $(OBJECTS)

.f90.o;
[TAB] $(COMP_FORTRAN) -c $.f90

.f.o:;
[Tab] $(COMP_FORTRAN) -c $.f
}}}

!!!2. Linking to external code or libraries
In this case, you will learn how to use pieces of code ("libraries") already written by others. As an example, we will use already existing external code to find the accuracy of the machine that you are using.

The most complicated part of the task is determining how to write your code so that it properly interfaces with the required library. All libraries have documentation, and you want to pay special attention to the arguments that are required and returned. In the case at hand:
* In C. Open with a text editor the file dlamch.c Line 54 describes the purpose of this library. Line 59 and following describes the arguments needed to interface with the library. Line 21 gives the entry point to the library, which indicates how you should call it from your code.
* In FORTRAN. Open with a text editor the file dlamch.f. The arguments are described in line 20 and following. How to call it can be determined from line 1.

Either:
* From github.umn.edu/jvinals/~PHYS4041-2018 download the subroutines/functions dlamch.f and lsame.f (FORTRAN), dlamch.c, lsame.c blaswrap.h (C) that are part of the LAPACK package. {{{git pull}}} should do it if you are set up correctly with github.
* Write the main code that uses this subroutines/functions to obtain information about the accuracy of the computer.
* To compile your code, e.g., dlam_drive.c or dlam_drive.f90, and produce the executable dlamch_drive type
** FORTRAN 
{{{
gfortran dlamch_drive.f90 -o dlamch_drive dlamch.f lsame.f
}}}
** C
{{{
gcc dlamch_drive.c -o dlamch_drive dlamch.c lsame.c /usr/lib/x86_64-linux-gnu/libf2c.so  -u MAIN__
}}}

Or:
* Alternatively, you can directly link to the LAPACK set of libraries that is already installed in your computer. The compiler will find and use the subroutines/functions that you need without any need to download them one by one. The entire LAPACK package (numerical linear algebra) already compiled is in the archive (CSELABS) /usr/lib/liblapack.so for FORTRAN and /usr/lib/liblapacke.so (C), and (Physics) /usr/lib64/liblapack.so (FORTRAN) and /usr/lib64/liblapacke.so (C). For example, if your main code were called dlamch_drive.f90, you would compile it on a Physics department system, for example, by typing:
{{{
gfortran dlamch_drive.f90 -o dlamch_drive /usr/lib64/liblapack.so
}}}
which is often abbreviated into
{{{
gfortran dlamch_drive.f90 -o dlamch_drive -L /usr/lib64/ -llapack
}}}
In C, you would type something similar. In the CSELABS systems, for example, you would type
{{{
gcc dlamch_drive.c  -o dlamch_drive -L /usr/lib -llapack
}}}
Create a makefile to automate the code build.
!!!!Combination of spline finding and ~Newton-Raphson method

In problem 2 of Homework three, you wrote a code that uses spline interpolation to find, given the data, the volume $v$ at which the interpolation gives $p=3.25$.

Recast the problem as finding the solution to $p(v) - 3.25 = 0$. Use the ~Newton-Raphson algorithm to accomplish this ''with'' the spline interpolant as the function.

Make sure you write the appropriate makefile to compile this set of codes.
!!!!1. Bessel function approximation
Bessel functions are trancendental functions that can be obtained by summing a series. However, there are recursion formulae that make the computation more efficient, although, as seen in class, recursion can introduce its own numerical problems. Use the following recursion formulae to write a code that computes the Bessel function $J_{l}$ with $l = 0, 1, 2, 3, 4$ and plot the results for $x \in (0,12)$. The recursion formulae ("upwards" or "downwards") are \begin{eqnarray}
J_{l+1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l-1}(x) \\
J_{l-1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l+1}(x) \\
J_{0}(x) = \frac{\sin x}{x} & \; & J_{1}(x) = \frac{\sin x - x \cos x}{x^{2}}
\end{eqnarray}  Using the following table to help evaluate the accuracy of your approximations,
|! x |! $J_{3}$ |! $J_{5}$ |! $J_{8}$ |
| 0.10 | 9. 518519719 10$^{-6}$ | 9.616310231 10$^{-10}$ | 2.901200102 10$^{-16}$ |
| 1.00 | 9.006581118 10$^{-3}$ | 9.256115862 10$^{-5}$ | 2.826498802 10$^{-8}$ |
| 10.0 | -3.949584498 10$^{-2}$ | -5.553451162 10$^{-1}$ | 1.255780236 10$^{-1}$ |
** Compute the Bessel functions required by using the upward recursion.
** Compare your solution to the table. Argue how errors might have propagated and what effect this may have on the accuracy of the solution.
** Chose an arbitrary starting function $J_{L}(x)$ with large $L$ and use the iteration downward. Compare the results to the upward iteration for the Bessel functions asked above. Note that the downward recursion relation is unchanged if you multiply both sides by an arbitrary function $f(x)$. As a consequence, for an arbitrary starting function, the sequence obtained is defined up to this factor. Let $J'_{0}(x)$ be the zero-th order obtained, and $J_{0}(x)$ the exact function as given. Compute the ratio $J_{0}(x)/J'_{0}(x)$ and scale ''all'' functions $J_{l}$ obtained by this ratio.
{{justifyright{
''Thursday, September 20, 2018''
}}}

A simple procedure to calculate numerically the //sin// function is through its Taylor series
$$
\sin(x) \approx \sin_{N}(x) = \sum_{l=0}^{l=N} (-1)^{l} \frac{x^{2l+1}}{(2l+1)!} = \sum_{l=0}^{N} A_{l}
$$
Write a code that that approximates $\sin(x)$ by $\sin_{13}(x)$ for 100,000 points uniformly distributed in $[ -\pi,\pi ]$, that is, by the first 14 terms in the Taylor series expansion.

As was the case in the homework, this is an alternating series, and it can be evaluated more efficiently through recursion. In particular, we have the result
$$
A_{l} = - \frac{x^{2}}{2l(2l+1)} A_{l-1} \quad {\rm and} \; A_{0} = x
$$
# Compute the same values of $sin(x)$ by this recursion, and compare them to what you obtained by summing the series directly (e.g., by plotting both).
# Compare the execution times of the two methods, and compare them to using directly the intrinsic function $sin(x)$ provided by your computer.
!!!!! Computing elapsed execution time in FORTRAN
{{{
CALL cpu_time(t0) ! Clock starts
! ... tasks to be timed ...
CALL cpu_time(t1) 
WRITE (6,*) 'Elapsed time (in seconds): ',t1-t0
}}}
!!!!! Computing elapsed execution time in C++
{{{
#include <time.h> 

clock_t time; // start clock time
double diff;  // interval time in sec
.....

time = clock();      // start the clock

// Include statements that you want to time

diff = (double)(difftime(clock(), time))/CLOCKS_PER_SEC; // elapsed time since start
}}}
!!Solution of a system of equations by using LAPACK
Lapack is an open source, widely used set of routines to accomplish a wide range of numerical linear algebra tasks. We will use it here to solve a linear system of the type $Ax =b$.
Information about Lapack:
* User Guide: http://www.netlib.org/lapack/lug/
* We will be using the routine dgesv. Detailed information at: http://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f.html#a5ee879032a8365897c3ba91e3dc8d512 This is the direct URL for this routine. A useful interface for specific routine usage is: http://www.netlib.org/lapack/explore-html/index.html
* C interface for Lapack: http://www.netlib.org/lapack/lapacke.html
!!!!Coding instructions - FORTRAN users
Computers in the Physics department (not CSELAB) have Lapack installed. You will have to ssh to a Physics host. Then when you compile, you need to link the Lapack libraries. All libraries can be linked at compilation time as
{{{
gfortran your_code.f90 -o your_code_executable -L /usr/lib64/ -llapack
}}}
!!!!Coding instructions - C users
Neither Physics computers nor CSELAB have the C version. We have uploaded the libraries into the class repository. The libraries are in liblapacke.a . Header files are also in the respository. Navigate to the directory in which you created the git repository for this class, and update it:
{{{
git pull
}}}
You would then link your code as
{{{
gcc/g++ your_code.c -o your_code_executable -llapack -lblas liblapacke.a -I include/
}}}
Note that the statement
{{{
#include "lapacke.h"
}}}
is needed in your code your_code.c. Note also that the function name is ~LAPACKE_dgesv, with the same calling arguments as the FORTRAN [[version|http://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f.html#a5ee879032a8365897c3ba91e3dc8d512]].
!!!!1. Solution of a simple system of linear equations
Write a code (and makefile) that using the Lapack routine dgesv solves the following system of equations:
\begin{eqnarray}
3 x_{1} + 2 x_{2} & = & 2 \\
2 x_{1} + 6 x_{2} & = & -8
\end{eqnarray}
!!!!2. Scaling analysis of the routine dgesv
We wish to examine how the execution time grows with $n$, the order of the system of equations to solve. In order to do so,
# Create a large matrix $A$ with random entries. CALL srand(int seed) (seed is an arbitrary integer) will initialize the random sequence, and rand() returns real random numbers uniformly distributed between 0 and 1. Create a matrix of size $(n,n)$ with random entries, and also a right hand side vector $b$. To make the problem better conditioned, add to the diagonal elemens of $A$ some positive quantity (e.g. 1).
# Solve the linear system of equations for increasing values of $n$. Use the time command to obtain the execution time. For example, if a.out is the executable, type instead 
{{{
time a.out
}}}
This will give you the execution time. Plot in a log log scale the execution time with $n$ and examine how the time scales with $n$ for large $n$.
Although simple, this method will also include the time involved in the other tasks in your code. For example, it will include the time spent setting up the matrix A. A more accurate timing can be obtained by using the function cpu_time (in fortran):
{{{
CALL cpu_time(t0)
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE(6,*) 'Elapsed time: ',t1-t0
}}}
In C you would use:
{{{
#include <time.h>
clock_t start, end;
double cpu_time_used;

start = clock();
… /* Do the work. */
end = clock();
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
}}}
We will use the LAPACK routine dgesv to solve a partial differential equation.

The heat diffusion equations in one dimension is,
$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}}$$
Assume we denote by $u_{i}^{n}=u(x_{i},t_{n})$ the solution at a discrete set of points in some interval of the $x$ axis, $x_{i} = i \Delta x$. $i = 0, \ldots N$, and corresponding to some discrete time $t_{n} = n \Delta t$, with $\Delta x$ and $\Delta t$ the discrete increments in space and time respectively. We are told that the differential equation can be written in these discretized variables as
$$ \frac{1}{\Delta t} \left( u_{i}^{n+1}-u_{i}^{n} \right) = \frac{D}{(\Delta x)^{2}} \left( u_{i+1}^{n+1} - 2 u_{i}^{n+1} + u_{i-1}^{n+1} \right) $$

By rearranging the terms of this equation, an iterative solution can be constructed that relates the solution at some time $t_{n}$, $u_{i}^{n}$, to the solution at the next time step $t_{n+1}$, $u_{i}^{n+1}$. The iteration is,
$$ -a u_{i-1}^{n+1} + (1+2a) u_{i}^{n+1} - a u_{i+1}^{n+1} = u_{i}^{n}$$
with $a = D \Delta t/(\Delta x)^{2}$. This iteration can be written in matrix form $Ax=b$ where $x=\{u_{i}^{n+1} \}$, $b=\{u_{i}^{n} \}$ and the matrix $A$ is $$A = \left( \begin{array}{cccccc} 1+2a & -a & 0 & \ldots& & 0 \\ -a & 1+2a & -a & 0 & \dots & 0 \\ 0 & -a & 1+2a & -a & \ldots & 0 \\ & \ldots & \ldots & \ldots & & \\ 0 & \ldots & 0 & & -a & 1+2a \end{array} \right)$$

As boundary conditions at the end points, we take fixed $u$ for all times: $u_{0} = 0 $ and $u_{N} = 0$.

Write a code that uses the LAPACK routine dgesv to solve the linear system of equations given which amounts to solving the partial differential equation iteratively: Solve for a number of time steps $n = 1, 2, \ldots$ with an initial condition at $t=0$ $u_{k}^{0}=1$ for $k = N/2$ and the remaining $u_{k}^{0} = 0$.

Plot the solution for a few time steps as it evolves in time.
{{justifyright{
''Thursday, September 27, 2018''
}}}

We will use the LAPACK routine dgesv to solve a partial differential equation.

The heat diffusion equations in one dimension is,
$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}}$$
Assume we denote by $u_{i}^{n}=u(x_{i},t_{n})$ the solution at a discrete set of points in some interval of the $x$ axis, $x_{i} = i \Delta x$. $i = 0, \ldots N$, and corresponding to some discrete time $t_{n} = n \Delta t$, with $\Delta x$ and $\Delta t$ the discrete increments in space and time respectively. We are told that the differential equation can be written in these discretized variables as
$$ \frac{1}{\Delta t} \left( u_{i}^{n+1}-u_{i}^{n} \right) = \frac{D}{(\Delta x)^{2}} \left( u_{i+1}^{n+1} - 2 u_{i}^{n+1} + u_{i-1}^{n+1} \right) $$

By rearranging the terms of this equation, an iterative solution can be constructed that relates the solution at some time $t_{n}$, $u_{i}^{n}$, to the solution at the next time step $t_{n+1}$, $u_{i}^{n+1}$. The iteration is,
$$ -a u_{i-1}^{n+1} + (1+2a) u_{i}^{n+1} - a u_{i+1}^{n+1} = u_{i}^{n}$$
with $a = D \Delta t/(\Delta x)^{2}$. This iteration can be written in matrix form $Ax=b$ where $x=\{u_{i}^{n+1} \}$, $b=\{u_{i}^{n} \}$ and the matrix $A$ is $$A = \left( \begin{array}{cccccc} 1+2a & -a & 0 & \ldots& & 0 \\ -a & 1+2a & -a & 0 & \dots & 0 \\ 0 & -a & 1+2a & -a & \ldots & 0 \\ & \ldots & \ldots & \ldots & & \\ 0 & \ldots & 0 & & -a & 1+2a \end{array} \right)$$

As boundary conditions at the end points, we take fixed $u$ for all times: $u_{0} = 0 $ and $u_{N} = 0$.

Write a code that uses the LAPACK routine dgesv to solve the linear system of equations given. This amounts to solving the partial differential equation iteratively over $n$. Solve for a number of time steps $n = 1, 2, \ldots$ with an initial condition at $t=0$ $u_{k}^{0}=1$ for $k = N/2$ and the remaining $u_{k}^{0} = 0$.

Plot the solution for a few time steps as it evolves in time.
The heat diffusion equations reads, in one dimension,
$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}}$$
Assume we denote by $u_{i}^{n}$ the solution at $x = i \Delta x$ corresponding to a time $t = n \Delta t$, with $\Delta x$ and $\Delta t$ the discrete increments in space and time. The differential equation is now discretized as
$$ \frac{1}{\Delta t} \left( u_{i}^{n+1}-u_{i}^{n} \right) = \frac{D}{(\Delta x)^{2}} \left( u_{i+1}^{n+1} - 2 u_{i}^{n+1} + u_{i-1}^{n+1} \right) $$

This iteration can be rewritten as:
$$ -a u_{i-1}^{n+1} + (1+2a) u_{i}^{n+1} - a u_{i+1}^{n+1} = u_{i}^{n}$$
with $a = D \Delta t/(\Delta x)^{2}$. We will also assume boundary conditions at the end points fixed for all times: $u_{0}$ and $u_{N}$.

Write a code that:
* Implements the iteration above in matrix form. Note that the matrix is tridiagonal.
* Implement the recursion relation given in class to solve the resulting tridiagonal system for one time step.
* Solve for a number of time steps in the case in which $u_{0} = u_{N} = 0$, and an initial condition $u_{k}^{0}=1$ for one arbitrary value of $k$ and the remaining $u^{0} = 0$.
!!! Solution of the biharmonic equation in two dimensions
We want to reuse the method that you have developed to solve the Poisson equation to solve the biharmonic equation instead
\begin{equation}
\nabla^{4} V(x,y) + b \nabla^{2} V(x,y) = -\rho(x,y) 
\end{equation}
where $b$ is a constant. By introducing the auxiliary function $\nabla^{2} V(x,y) = g(x,y)$, we decompose the original problem into two problems. First, the biharmonic equation can be written as $\nabla^{2}g+bg = -\rho$, which can solved for $g(x,y)$. Then, the definition of $g(x,y)$ is a Poisson equation determining $V(x,y)$.

Use the same discretization scheme and iterative algorithm that you have used in the homework for the Poisson equation (with $V_{ij} = V(x_{i},y_{j})$, $x_{i} = i\Delta x, y_{j} =j \Delta y$):
\begin{equation}
\frac{1}{\Delta x^{2}} \left( V_{i+1,j} - 2 V_{i,j} + V_{i-1,j} + V_{i,j+1} - 2 V_{i,j} + V_{i,j-1} \right) = - \rho_{i,j} \quad i = 1, ..., N; ~ j = 1, ... , M
\end{equation}
# First, modify the algorithm slightly as the equation for $g(x,y)$ is not the Poisson but rather the Helmholtz equation. Use the same iterative method as in the homework to solve this equation.
# Once $g(x,y)$ has been obtained, solve the Poisson equation with the same method to find $V(x,y)$.
{{justifyright{
''Thursday, October 4, 2018''
}}}
!!! Shear stress in turbulent flow
Theory suggests that for turbulent flow in a circular tube of radius $R=1$ the time smoothed velocity field $v_{z}$ along the tube axis is
$$v_{z} = v_{max} \left( 1 - \frac{r}{R} \right)^{\alpha}$$ 
where $r$ is the radial coordinate inside the tube, and $\alpha$ an unknown exponent. The following measurements have been performed on the velocity as a function of $r$
| $r$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 |
| $v_{z}$ | 1.13 | 1.10 | 1.05 | 0.90 | 0.82 | 0.65 | 0.43 | 0.25 | 0.12 |

We are interested in computing from the data the shear stress in the fluid $T$, which in cylindrical coordinates is given by
$$ T = \left( \frac{d^{2} v_{z}}{dr^{2}} + \frac{1}{r} \frac{d v_{z}}{dr} \right)$$

Write a code that computes numerically the shear stress by
# Directly approximating the derivatives by finite differences of the data.
# Interpolate a cubic spline to the data to find $T$ in the interval $r \in (0.1,0.9)$
!!!!1. Continuation of Lab Session 5
* Finish writing the code to solve the diffusion equation as suggested in [[LabSession 5]].
* Integrate que equation forward in time for the following initial conditions:
** A rectangular temperature peak at the center of the domain.
** Two rectangular temperature peaks somewhere in the domain.
** Fixed boundary contitions $u_{0}=2$ and $u_{N} = 1$, with all other $u = 0$ at $t=0$.
** In all three cases, plot the resulting solutions for several times.

!!!!2. Integrals over unbounded domains
Time permitting, try to evaluate the following integral over an unbounded domain,
$$ I = \int_{\pi}^{\infty} (s + x)^{-1/3} \sin x dx$$ by using the trapezoidal rule.
@@color(blue):Hint:@@ Compute $I_{n}^{(0)} = \int_{\pi}^{n \pi} (s + x)^{-1/3} \sin x dx$ by using the trapezoidal rule for a range of values of $n$. The resulting sequence $I_{n}^{(0)}$ is alternating and converges very slowly with increasing $n$. To accelerate convergence define $I_{n}^{(1)} = (1/2)(I_{n}^{(0)} + I_{n-1}^{(0)})$. Write a code that iterates this averaging $I_{n}^{(k+1)} = (1/2)(I_{n}^{(k)} + I_{n-1}^{(k)})$ to estimate the value of the integral.
!!! Integrals over unbounded domains
Evaluate the following integral over an unbounded domain,
$$ I = \int_{\pi}^{\infty} (s + x)^{-1/3} \sin x dx$$ by using the trapezoidal rule.
@@color(blue):Hint:@@ Compute $I_{n}^{(0)} = \int_{\pi}^{n \pi} (s + x)^{-1/3} \sin x dx$ by using the trapezoidal rule for a range of values of $n$. The resulting sequence $I_{n}^{(0)}$ is alternating and converges very slowly with increasing $n$. To accelerate convergence define $I_{n}^{(1)} = (1/2)(I_{n}^{(0)} + I_{n-1}^{(0)})$. Write a code that iterates this averaging $I_{n}^{(k+1)} = (1/2)(I_{n}^{(k)} + I_{n-1}^{(k)})$ to estimate the value of the integral.

!!! Eigenvalue problems for differential equations
Write a code that computes the smallest eigenvalue $\lambda$ for which the problem $$
\frac{d}{dx} \left[ (1+x^{2}) \frac{dy}{dx} \right] + \lambda y = 0, \quad\quad y(-1)=y(1)=0
$$
has a nontrivial solution.

In order to do it, use the following discretization scheme, $$ 
\frac{d}{dx} \left[ p(x) \frac{dy}{dx} \right] = \frac{1}{h} \left[ p_{n+1/2}\left( \frac{y_{n+1}-y_{n}}{h} \right) - p_{n-1/2} \left( \frac{y_{n}-y_{n-1}}{h} \right) \right] + c_{1}(x) h^{2} + c_{2} h^{4} + \ldots ,
$$
where we have also included the asymptotic error of the discretization. Once the differential equation is discretized, find the smallest eigenvalue of the resulting matrix.
{{justifyright{
''Thursday, October 11, 2018''
}}}

We investigate here a common technique of code development: the use of preprocessors. They "pre process" the source code @@color(red):before@@ compilation. By introducing special commands ("directives") the preprocessor can act on the code prior to compilation,  and hence modify the actual code send to the compiler. For that reason, they add a great deal of flexibility in code management and reduce errors by simplifying version maintenance.

!!!Preprocessing (cpp)
As an example, we will write a code that integrates an ordinary differential equation using either forward or backward Euler methods, and we will use a pre-processor to determine which algorithm is used. For C, C++, and FORTRAN code we will use the C preprocessor: cpp (man cpp for help).

cpp processes special statements called directives. They are lines in the source code (in C, C++, and FORTRAN) that begin with #. For example
{{{
# define PI 3.141592
}}}
would literally replace every instance of PI (case sensitive) in a given source code by 3.141592. As an example in FORTRAN, create the code test.f90:
{{{
PROGRAM test
# define PI 3.141592

x = PI
write(6,*) x

END PROGRAM
}}}
In C++ you would use:
{{{
#include <iostream>
#define PI 3.141592

using namespace::std;

int main(){
  cout << "pi = " << PI << endl;
  return 0;
}
}}}

To run the code through the preprocessor, you would type
{{{
cpp -P test.f90
}}}
The standard ouput will give you the "preprocessed" code. Note how PI has been replaced by its value. To store the preprocessed code in a different file (e.g., test_processed.f90), you would use
{{{
cpp -P -o test_processed.f90 test.f90
}}}
This new file does not contain PI, rather its numerical value. This is a useful feature, for example, in allocating array dimensions without having to edit the source code (even more useful if the modification would involve a large number of modules and subroutines). More about this later.

Another useful feature of defined variables is that they can be used to control code. For example:
{{{
# define LINEAR
....
# if defined(LINEAR)
Execute some instructions
# endif

# if defined(NONLINEAR)
Some other instructions
# endif
}}}
The instructions in the first if block will only appear in the preprocessed code if LINEAR is defined above, but not otherwise. The first line (#define LINEAR) can in effect be used to control whether the LINEAR or NONLINEAR part of the code will be compiled. Try an example yourself, and type cpp -w test.f90 to see the effect.

A third useful directive is
{{{
# include "config.h"
}}}
which includes the file config.h before the preprocessig is done. For example, one could have a file config.h containing:
{{{
# define NSIZE 200
# define LINEAR
}}}
with the code test.f90 including now:
{{{
PROGRAM test

# include "config.h"

REAL(8) x(NSIZE)

# if defined (LINEAR)
do i=1,NSIZE
x(i)=i
END DO
# endif

# if defined(NONLINEAR)
do i=i,NSIZE
x(i)=i**2
END DO
# endif

END PROGRAM
}}}

And the corresponding makefile could read something like
{{{
test: test.f90 config.h
	cpp -P -o test_cpp.f90 test.f90
	gfortran -o test test_cpp.f90
	rm test_cpp.f90
}}}

@@color(red):NOTE:@@ In practice, compilers know how to interpret pre-processor directives without the need to create intermediate files (the files test_cpp.f90 above). In C/C++ the preprocessor is used automatically, and there is nothing for you to do anything at compilation time. In FORTRAN you can either use gfortran -cpp test.f90 or gfortran test.F90. In the first case, the compiler is alerted via -cpp to run the code with the preprocessor. In the second, the convention is used that if the extension has a capital F (test.F90) not a lower case (test.f90), the compiler will use the preprocessor first.

Use this method to write a code that solves the ordinary differential equation
$$ \frac{dy}{dt} = \frac{\sin y}{\sqrt{y}} $$
with initial condition at $t=0$ $y=1$ up to $t=1$. Write a code that can use either the forward or backward Euler method and use the pre-processor to choose which method to use. For example, create a config.h file that contains # define FORWARD or # define BACKWARD and the corresponding sections of the main code (e.g. euler.f90 or euler.c) to use one algorithm or the other depending on the contents of the file config.h 
!!!!Boundary value problem
We wish to solve the following second order equation in one dimension
$$y^{\prime\prime} = -\frac{4}{x} y^{\prime} + \frac{2}{x^{2}} y - \frac{2}{x^{2}} \ln x$$
in the interval $1 \le x \le 2$, and with boundary conditions $y(1) = -1/2$, $y(2) = \ln 2$.
# Write the second order equation as a system of two first order equations, each with one boundary condition.
# Use a method of your choice to discretize and solve the system of two equations as if it were an initial value problem.
# We will use a "shooting method" to solve the boundary value problem. Set up the solution method by assuming one missing boundary condition at $x = 1$ and integrate the system up to $x=2$.
# Set up an interative method so that changes in the assumed boundary condition at $x=1$ lead, after convergence, to the correct boundary condition at $x=2$, and hence to the solution of the boundary value problem.

A simple minded solution. The code bvp.f90:
{{{
PROGRAM bvp

p0_old=3.0
dp=0.01
CALL euler(0.01,-0.5,p0_old,yn1,pn1)
p0_plus=p0_old+dp
CALL euler(0.01,-0.5,p0_plus,yn2,pn2)

p0_new=p0_old-(yn1-log(2.0))*dp/(yn2-yn1)

CALL euler(0.01,-0.5,p0_new,yn3,pn3)

WRITE(6,*) p0_old,yn1-log(2.0),p0_new,yn3-log(2.0)

END PROGRAM
}}}
whis uses the subroutine euler.f90
{{{
SUBROUTINE euler(dx,y0,p0,yn,pn)

open(1,file='output.dat')

yn=y0
pn=p0

niter=(2.0-1.0)/dx

do i=1,niter
ysave=yn
x=1.0+i*dx
yn=yn+dx*pn
pn=pn+dx*(-4.0*pn/x+2.0*ysave/x**2-2.0*log(x)/x**2)
write(1,*) x,yn,pn
end do

CLOSE(1)

RETURN
END SUBROUTINE
}}}
and the makefile to put it all together
{{{
bvp: bvp.f90 euler.o
        gfortran -o bvp bvp.f90 euler.o

euler.o: euler.f90
        gfortran -c euler.f90
}}}
Numerical integration of a predator-prey model. Consider two species of fish that coexist in the same lake. Let $x$ be the 
!!! Numerical integration of a predator-prey model. 

Consider two species of fish that coexist in the same lake. Let $x$ be  number of small fish of type A which only eat grass from the bottom, and $y$ the number of large fish of type B which only eat small fish of type A, but no grass. Alone, the population of A might be expected to grow exponentially, whereas B alone (no A to eat) would die exponentially. However, the two species interact: whenever a fish B encounters a fish A, it will eat it. This interaction contributes to the increase in the population of B and a decrease in A, both rates being proportional to the quantity $xy$. A general mathematical model of this system is $$ \frac{dx}{dt} = \epsilon_{A} x - \gamma_{A} xy $$
$$ \frac{dy}{dt} = - \epsilon_{B} y + \gamma_{B} xy $$  where $\epsilon_{A}$ is the growth rate of A alone, $\epsilon_{B}$ the extinction of B alone, and $\gamma_{A}$ and $\gamma_{B}$ the interaction coefficients. All four are assumed to be constant. Write a code that integrates this system of differential equations with $\epsilon_{A} = \epsilon_{B} = \gamma_{A} = \gamma_{B} =1$, and initial conditions at $t=0$ $x=0.1$ and $y=0.1$. 
## Integrate the system up to $t=10$.
## Plot $y$ versus $x$ and provide an interpretation of the resulting graph.
{{justifyright{
''Thursday, October 18, 2018''
}}}
!!! Boundary value problem
We wish to solve the following second order equation in one dimension
$$y^{\prime\prime} = -\frac{4}{x} y^{\prime} + \frac{2}{x^{2}} y - \frac{2}{x^{2}} \ln x$$
in the interval $1 \le x \le 2$, and with boundary conditions $y(1) = -1/2$, $y(2) = \ln 2$.
# Write the second order equation as a system of two first order equations, each with one boundary condition.
# Use a method of your choice to discretize and solve the system of two equations as if it were an initial value problem.
# We will use a "shooting method" to solve the boundary value problem. Set up the solution method by assuming one missing boundary condition at $x = 1$ and integrate the system up to $x=2$.
# Set up an interative method so that changes in the assumed boundary condition at $x=1$ lead, after convergence, to the correct boundary condition at $x=2$, and hence to the solution of the boundary value problem.
We investigate here a common technique of code development: the use of preprocessors. They "pre process" the source code @@color(red):before@@ compilation. By introducing special commands ("directives") the preprocessor can act on the code prior to compilation,  and hence modify the actual code send to the compiler. For that reason, they add a great deal of flexibility in code management and reduce errors by simplifying version maintenance.

!!!Preprocessing (cpp)
As an example, we will write a code that integrates an ordinary differential equation using either forward or backward Euler methods, and we will use a pre-processor to determine which algorithm is used. For both C and FORTRAN code we will use the C preprocessor: cpp (man cpp for help).

cpp processes special statements called directives. They are lines in the source code (both C and FORTRAN) that begin with #. For example
{{{
# define PI 3.141592
}}}
would literally replace every instance of PI (case sensitive) in a given source code by 3.141592. As an example, create the code test.f90:
{{{
PROGRAM test
# define PI 3.141592

x = PI
write(6,*) x

END PROGRAM
}}}
To run the code through the preprocessor, you would type
{{{
cpp -w test.f90
}}}
The standard ouput will give you the "preprocessed" code. Note how PI has been replaced by its value. To store the preprocessed code in a different file (e.g., test_processed.f90), you would use
{{{
cpp -w -o test_processed.f90 test.f90
}}}
This new file does not contain PI, rather its numerical value. This is a useful feature, for example, in allocating array dimensions without having to edit the source code (even more useful if the modification would involve a large number of modules and subroutines). More about this later.

Another useful feature of defined variables is that they can be used to control code. For example:
{{{
# define LINEAR
....
# if defined(LINEAR)
Execute some instructions
# endif

# if defined(NONLINEAR)
Some other instructions
# endif
}}}
The instructions in the first if block will only appear in the preprocessed code if LINEAR is defined above, but not otherwise. The first line (#define LINEAR) can in effect be used to control whether the LINEAR or NONLINEAR part of the code will be compiled. Try an example yourself, and type cpp -w test.f90 to see the effect.

A third useful directive is
{{{
# include "config.h"
}}}
which includes the file config.h before the preprocessig is done. For example, one could have a file config.h containing:
{{{
# define NSIZE 200
# define LINEAR
}}}
with the code test.f90 including now:
{{{
PROGRAM test

# include "config.h"

REAL(8) x(NSIZE)

# if defined (LINEAR)
do i=1,NSIZE
x(i)=i
END DO
# endif

# if defined(NONLINEAR)
do i=i,NSIZE
x(i)=i**2
END DO
# endif

END PROGRAM
}}}

And the corresponding makefile could read something like
{{{
test: test.f90 config.h
	cpp -w -o test_cpp.f90 test.f90
	gfortran -o test test_cpp.f90
	rm test_cpp.f90
}}}

Use this method to write a code that solves the ordinary differential equation
$$ \frac{dy}{dt} = \frac{\sin y}{\sqrt{y}} $$
with initial condition at $t=0$ $y=1$ up to $t=1$. Write a code that can use either the forward or backward Euler method and use the pre-processor to choose which method to use. For example, create a config.h file that contains # define FORWARD or # define BACKWARD and the corresponding sections of the main code (e.g. euler.f90 or euler.c) to use one algorithm or the other depending on the contents of the file config.h 

!!!!Boundary value problem
We wish to solve the following second order equation in one dimension
$$y^{\prime\prime} = -\frac{4}{x} y^{\prime} + \frac{2}{x^{2}} y - \frac{2}{x^{2}} \ln x$$
in the interval $1 \le x \le 2$, and with boundary conditions $y(1) = -1/2$, $y(2) = \ln 2$.
# Write the second order equation as a system of two first order equations, each with one boundary condition.
# Use a method of your choice to discretize and solve the system of two equations as if it were an initial value problem.
# We will use a "shooting method" to solve the boundary value problem. Set up the solution method by assuming one missing boundary condition at $x = 1$ and integrate the system up to $x=2$.
# Set up an interative method so that changes in the assumed boundary condition at $x=1$ lead, after convergence, to the correct boundary condition at $x=2$, and hence to the solution of the boundary value problem.

A simple minded solution. The code bvp.f90:
{{{
PROGRAM bvp

p0_old=3.0
dp=0.01
CALL euler(0.01,-0.5,p0_old,yn1,pn1)
p0_plus=p0_old+dp
CALL euler(0.01,-0.5,p0_plus,yn2,pn2)

p0_new=p0_old-(yn1-log(2.0))*dp/(yn2-yn1)

CALL euler(0.01,-0.5,p0_new,yn3,pn3)

WRITE(6,*) p0_old,yn1-log(2.0),p0_new,yn3-log(2.0)

END PROGRAM
}}}
whis uses the subroutine euler.f90
{{{
SUBROUTINE euler(dx,y0,p0,yn,pn)

open(1,file='output.dat')

yn=y0
pn=p0

niter=(2.0-1.0)/dx

do i=1,niter
ysave=yn
x=1.0+i*dx
yn=yn+dx*pn
pn=pn+dx*(-4.0*pn/x+2.0*ysave/x**2-2.0*log(x)/x**2)
write(1,*) x,yn,pn
end do

CLOSE(1)

RETURN
END SUBROUTINE
}}}
and the makefile to put it all together
{{{
bvp: bvp.f90 euler.o
        gfortran -o bvp bvp.f90 euler.o

euler.o: euler.f90
        gfortran -c euler.f90
}}}
{{justifyright{
''Thursday, October 25, 2018''
}}}
!!! ~Runge-Kutta integration of a system of ordinary differential equations
We wish to calculate the vertical density distribution of a self-gravitating plane distribution of gas in isothermal, hydrostatic equilibrium by using the ~Runge-Kutta solver of order $4^{th}$. The model is defined by by the following two, coupled ordinary differential equations,
$$
\frac{d \rho(z)}{dz} = \frac{g}{c^{2}} \rho ,
$$
$$
\frac{d g(z)}{dz} = - 4 \pi G \rho ,
$$
where $\rho(z)$ is the mass density and $g(z)$ the self-consistently determined intensity of the gravitational field. $c^{2}$ and $G$ are known constants. For convenience, you can set $c^{2} = 2 \pi G = 1$.

# Write an algorithm using the ~Runge-Kutta method of order 4 to solve this system of equations with boundary conditions $\rho(z=0) = \rho_{0} = 1$ and $g(z=0) =0$ up to $z_{max} = 10$. A discretization step $\Delta z = 0.4$ should be sufficient.
# For reference, compare your solution to the known analytic solution given by, $$
\rho(z) = \rho_{0} \; {\rm sech}^{2}(z/H) \quad g(z) = -4 \pi G \Sigma \; {\rm tanh}(z/H)
$$
where $H = c^{2}/(2 \pi G \Sigma)$ and $\Sigma = \int_{0}^{\infty} \rho(z) dz$. You can estimate $\Sigma$ by using the trapezoidal rule on the numerically determined density profile.
Refresh your git repository with updated tridiagonal routines. Go to the directory in which you have the respository and type
{{{
git pull
}}}

!!! Alternating Direction Implicit solution of Poisson's equation in two dimensions
As we have seen in class, a simple discretization of Poisson's equations on a two dimensional square grid of spacing $\Delta x$ is
$$u_{i+1 j} + u_{i-1 j} + u_{i j-1} + u_{i j+1} - 4 u_{ij} = - \Delta x^{2} \rho_{ij}$$
In order to develop and ADI algorithm, we split it as follows
$$\left( u_{i+1 j} + u_{i-1 j} - 2 u_{ij} \right) + \left( u_{i j-1} + u_{i j+1} - 2 u_{ij} \right) = - \Delta x^{2} \rho_{ij}$$
In the first parenthesis the index $j$ is constant, whereas in the second parenthesis the index $i$ is constant.

!!!!First half step
We define the first half step of the ADI method assuming implicit updating in $i$ only:
$$u_{i+1 j}^{n+1/2} + u_{i-1 j}^{n+1/2} - 2 u_{ij}^{n+1/2} = - \left( u_{i j-1}^{n} + u_{i j+1}^{n} - 2 u_{ij}^{n}  \right) - \Delta x^{2} \rho_{ij}$$

@@color(maroon):If you now think that $j$ is fixed@@, then for each $j$ this equations is just a standard tri-diagonal system - or in effect $n$ tridiagonal systems, one for each $j$. The right hand side depends on $\rho_{ij}$ but also on the previous solution $u_{ij}$. Note that the right hand side involves different values of $j$, but since they are at an earlier time, they are known. 

We also need to specify boundary conditions. For now, let us take $u = 0$ on all boundaries. Then, in the notation given in class, the elements of the tridiagonal matrix are:
\begin{eqnarray}
b(1) = b(n) = 1 & \quad & b(2) ... b(n-1) = -2 \\
a(n) = 0 & \quad & a(2) ... a(n-1) = 1\\
c(1) = 0 & \quad & c(2) ... c(n-1) = 1
\end{eqnarray} 

Write a loop over $j$ so that for each $j$ the tridiagonal system above is solved, leading to $u_{ij}^{n+1/2}$. Use the routine provided: tridag.f or tridag.c.

!!!!Second half step
We now reverse the implicit direction for the second half step:
$$u_{i j+1}^{n+1} + u_{i j-1}^{n+1} - 2 u_{ij}^{n+1} = - \left( u_{i+1 j}^{n+1/2} + u_{i-1 j}^{n+1/2} - 2 u_{ij}^{n+1/2}  \right) - \Delta x^{2} \rho_{ij}$$
@@color(maroon):If you now think of $i$ being fixed instead@@, for each $i$ we have a tridiagonal system on $j$. So in effect, we have again $n$ tridiagonal systems. Furthermore, note that the form of the tridiagonal system is identically the same as above.

Write a second loop to solve all $n$ tridiagonal systems. 

This completes one iteration in the solution to Poisson's equation for any given $\rho_{ij}$ and boundary conditions. The process can be iterated until convergence. Sample code:
{{{
PROGRAM poisson

PARAMETER (n=100)

REAL(8) :: u(n,n),un(n,n),rho(n,n),r(n,n)
REAL(8) :: a(n),b(n),c(n)

dx=0.1

do i=1,n
do j=1,n
rho(i,j)=-dx**2*sqrt(FLOAT(i)**2+FLOAT(j)**2)
u(i,j)=0.0
un(i,j)=0.0
end do
end do

a(n)=0.0
do i=1,n-1
a(i)=1.0
end do

b(1)=1.0
b(n)=1.0
do i=2,n-1
b(i)=-2.0
end do

c(1)=0.0
do i=2,n-1
c(i)=1.0
end do

DO j=2,n-1
DO i=1,n
r(i,j)=rho(i,j)-(u(i,j+1)+u(i,j-1)-2.0*u(i,j))
END DO
call tridag(a,b,c,r(:,j),un(:,j),n)
END DO

DO i=2,n-1
DO j=1,n
r(i,j)=rho(i,j)-(un(i+1,j)+un(i-1,j)-2.0*un(i,j))
END DO
call tridag(a,b,c,r(i,:),u(i,:),n)
END DO

WRITE(6,*) 'One iteration completed'

END PROGRAM
}}}
!!! Alternating Direction Implicit (ADI) method for convection-diffusion equation
The temperature distribution in a porous rock in the presence of ground water flow directed along the $x$ axis is described by the equation
$$ v \frac{\partial T}{\partial x} = \kappa \left( \frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} \right) $$
where $v$ is the constant water velocity and $\kappa$ is the thermal diffusivity coefficient.

Each iteration of the ADI method involves two consecutive steps, each being semi-implicit in alternating directions. In this case, considering central differences for the spatial derivatives, the two steps are $$
(1 -\beta) T_{i+1,j}^{n+1/2} + 2 \beta T_{i,j}^{n+1/2} - (1+\beta)  T_{i-1,j}^{n+1/2} = \beta \left(   T_{i,j+1}^{n}  -2 T_{i,j}^{n}  + T_{i,j-1}^{n} \right) 
$$  
and $$
\beta \left(  T_{i,j+1}^{n+1}  -2 T_{i,j}^{n+1}  + T_{i,j-1}^{n+1} \right) = (1-\beta)   T_{i+1,j}^{n+1/2} + 2   \beta T_{i,j}^{n+1/2} - (1+\beta)  T_{i-1,j}^{n+1/2}
$$ 
with $\beta = 2 \kappa/v \Delta x$ ($\Delta y = \Delta x$ here). The first step runs over all values of $j = 1, \ldots n$, and the second step runs over all values of $i = 1, \ldots, n$. Both are a tri-diagonal system of linear equations. 

Sketch the solution to the two systems of linear equations defined by these two iterations by using the routine tridag (FORTRAN or C) uploaded to the class repository. As you will see, the routine takes as input three arrays:  $b_{k}, k=1, \ldots N = n \times n \;$ contains the diagonal elements in the matrix representing the system, $a_{k}, k =2, \ldots, N$ the line immediately below the diagonal, and $c_{k}, k=1, \ldots N-1$ the line immediately above the diagonal. The right hand side is the array $r_{k}, k=1, \ldots N$.

In order to do so, map the two ADI steps given above to the corresponding arrays $a,b,c, r$. Then, given the output $u_{k}, k=1, \ldots N$ of tridag, map it back to the temperature $T_{ij}$.
{{justifyright{
''Thursday, November 1, 2018''
}}}
!!! Soliton solution
We wish to find nonlinear solitary wave solutions ("solitons") of the following equation $$
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \delta^{2} \frac{\partial^{3} u}{\partial x^{3}} = 0.
$$
where we choose $\delta = 0.022$ and take as initial condition $u(x,0) = \cos (\pi x)$, $0 \le x \le 2$. We will also impose boundary conditions that $u$ is periodic in $[0,2]$ for all time.

Develop a code that uses the ~Lax-Wendroff method to integrate the equation in time, and examine the emergence of the soliton solution. Note that you first have to write the equation given in conservative form $\partial_{t} u + \partial_{x}J(u) = 0$. Integrate the equation up to a maximum time $t=3.6/\pi$.
/***
|''Name:''|LoadRemoteFileThroughProxy (previous LoadRemoteFileHijack)|
|''Description:''|When the TiddlyWiki file is located on the web (view over http) the content of [[SiteProxy]] tiddler is added in front of the file url. If [[SiteProxy]] does not exist "/proxy/" is added. |
|''Version:''|1.1.0|
|''Date:''|mar 17, 2007|
|''Source:''|http://tiddlywiki.bidix.info/#LoadRemoteFileHijack|
|''Author:''|BidiX (BidiX (at) bidix (dot) info)|
|''License:''|[[BSD open source license|http://tiddlywiki.bidix.info/#%5B%5BBSD%20open%20source%20license%5D%5D ]]|
|''~CoreVersion:''|2.2.0|
***/
//{{{
version.extensions.LoadRemoteFileThroughProxy = {
 major: 1, minor: 1, revision: 0, 
 date: new Date("mar 17, 2007"), 
 source: "http://tiddlywiki.bidix.info/#LoadRemoteFileThroughProxy"};

if (!window.bidix) window.bidix = {}; // bidix namespace
if (!bidix.core) bidix.core = {};

bidix.core.loadRemoteFile = loadRemoteFile;
loadRemoteFile = function(url,callback,params)
{
 if ((document.location.toString().substr(0,4) == "http") && (url.substr(0,4) == "http")){ 
 url = store.getTiddlerText("SiteProxy", "/proxy/") + url;
 }
 return bidix.core.loadRemoteFile(url,callback,params);
}
//}}}
[[Home]]
[[Syllabus]]
[[Class Notes]]
[[Homeworks and Labs]]
[[Additional Reading]]
[[Tag Cloud]]

[[WelcomeToTiddlyspot]] [[GettingStarted]]
{{justifyright{
''Due Wednesday, October 22, 2014''
}}}

!!!!Shear stress in turbulent flow
Theory suggests that for turbulent flow in a circular tube of radius $R=1$ the time smoothed velocity field $v_{z}$ along the tube axis is
$$v_{z} = v_{max} \left( 1 - \frac{r}{R} \right)^{\alpha}$$ 
where $r$ is the radial coordinate inside the tube, and $\alpha$ an unknown exponent. The following measurements have been performed on the velocity as a function of $r$
| $r$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 |
| $v_{z}$ | 1.13 | 1.10 | 1.05 | 0.90 | 0.82 | 0.65 | 0.43 | 0.25 | 0.12 |

We are interested in computing from the data the shear stress in the fluid $T$, which in cylindrical coordinates is given by
$$ T = \left( \frac{d^{2} v_{z}}{dr^{2}} + \frac{1}{r} \frac{d v_{z}}{dr} \right)$$

Write a code that computes numerically the shear stress by
# Directly approximating the derivatives by finite differences of the data.
# Interpolate a cubic spline to the data to find $T$ in the interval $r \in (0.1,0.9)$

!!!!Differential equation solution
We wish to solve the equation
$$ \frac{dy}{dt} = \frac{\sin y}{\sqrt{y}} $$
with initial condition at $t=0$ $y=1$, with a predictor-corrector method. For the explicit part we will use the forward Euler method, whereas the implicit part of the algorithm uses the trapezoidal rule
$$ y_{i+1} = y_{i} + \frac{\Delta t}{2} \left[ f(t_{i},y_{i}) + f(t_{i+1},y_{i+1}) \right] $$
where $f(t,y)$ is the right hand side of the differential equation.
# Take $\Delta t = 0.1$ and write a code that calculates the solution up to $t=1$. For the implicit part, set up a fixed point iteration at each step to find the solution up to six decimal places.
# Derive the convergence condition for the fixed point iteration used, and give an upper bound for the time step $\Delta t$ that is necessary for the iteration to converge. 

[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/PHYS4041_MidTerm.pdf]]
{{justifyright{
''Due Friday, October 21, 2016, 5:00 p.m.''
}}}
We need to find the solution of the following integral equation in the interval $(0,1)$,$$
\int_{0}^{1} dx e^{-(x-y)^{2}/\epsilon^{2}} \frac{df(x)}{dx} = 2 \epsilon \sqrt{\pi} ( y - 3 y^{2} + 2 y^{3}),
$$
where the function $f(x)$ needs to satisfy the boundary conditions $f(0)=f(1)=0$. Also, $\epsilon = 0.1$.

Discretize $x_{j}= j h, j = 0, \ldots n$ for some $h$, so that $x_{0} = 0$ and $x_{n} = 1$, and $y_{i} = i h, i = 0, \ldots n$. Consistent with this discretization, one defines $f_{j} = f(x_{j})$. The boundary conditions are then $f_{0} = f_{n} = 0$.

# Write the discretized version of the integral equation by using the trapezoidal rule for the integral and central differences for the derivative at the inner points, the forward difference at $x=0$ and the backward difference at $x=1$.
# Use the discretized version to write the code that will find the solution $f(x)$ in $(0,1)$.
{{justifyright{
''Due Friday, October 19, 2018, 5:00 p.m.''
}}}

@@color(red):''NOTE:''@@ Just to re iterate that you are expected to work //independently// when solving the mid term exam. You can submit the exam electronically as you have done for the homework.

!!! Question 1. Numerical integration over an unbounded domain (50 points).
The integration methods discussed in class consider definite integrals over a finite domain $(a,b)$. However, it is possible to either transform the integrands or to analyze convergence with the limits of integration to evaluate integrals over unbounded domains. Propose a suitable plan, and write the associate code to evaluate the following integrals for $s = 0.5, 0.6, 0.7, ..., 3.0$. 
## $$\int_{0}^{\infty} \left( x^{3} + sx \right)^{-1/2} dx,$$  
## $$\int_{0}^{\infty} \left( x^{2} + 1 \right) ^{-1/2} e^{-sx} dx$$ 

!!! Question 2. Two predators, one prey model (50 points).
Consider two predator species of numbers $P_{1}$ and $P_{2}$ that prey on the same species of number $p$. A possible model of evolution in time $t$ given the interaction between the three species is,
$$
\frac{dp}{dt} = ap \left( 1 - \frac{p}{K} \right) - \left( b_{1}P_{1} + b_{2}P_{2} \right) p,
$$
$$
\frac{dP_{1}}{dt} = \epsilon_{1} b_{1} pP_{1} - m_{1} P_{1}, \quad\quad \frac{dP_{2}}{dt} = \epsilon_{2} b_{2} pP_{2} - m_{2} P_{2}.
$$
Consider the following values of the model parameters: $a = 0.2, K = 1.7, b_{1} = 0.1, b_{2} = 0.2, m_{1} = m_{2} = 0.1, \epsilon_{1}=1.0, \epsilon_{2} = 2.0$.
# Integrate the model forward in time to determine the populations as a function of time, given the initial conditions $p(0) = P_{2}(0) = 1.7$ and $P_{1}(0) = 1.0$.
# To understand the time dependence of the model, it may help to not only plot the functions $p(t), P_{1}(t)$ and $P_{2}(t)$ as a function of time, but also against each other.
# Based on the results obtained for the evolution, can two predators that share the same prey coexist ? You may want to change the parameters of the second predator and examine the changes induced in the evolution.
# Search for chaotic (random) behavior by varying the growth rates $a, \epsilon_{1}$ and $\epsilon_{2}$.

[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2018/PHYS4041_MidTerm.pdf]]
/***
|''Name:''|PasswordOptionPlugin|
|''Description:''|Extends TiddlyWiki options with non encrypted password option.|
|''Version:''|1.0.2|
|''Date:''|Apr 19, 2007|
|''Source:''|http://tiddlywiki.bidix.info/#PasswordOptionPlugin|
|''Author:''|BidiX (BidiX (at) bidix (dot) info)|
|''License:''|[[BSD open source license|http://tiddlywiki.bidix.info/#%5B%5BBSD%20open%20source%20license%5D%5D ]]|
|''~CoreVersion:''|2.2.0 (Beta 5)|
***/
//{{{
version.extensions.PasswordOptionPlugin = {
	major: 1, minor: 0, revision: 2, 
	date: new Date("Apr 19, 2007"),
	source: 'http://tiddlywiki.bidix.info/#PasswordOptionPlugin',
	author: 'BidiX (BidiX (at) bidix (dot) info',
	license: '[[BSD open source license|http://tiddlywiki.bidix.info/#%5B%5BBSD%20open%20source%20license%5D%5D]]',
	coreVersion: '2.2.0 (Beta 5)'
};

config.macros.option.passwordCheckboxLabel = "Save this password on this computer";
config.macros.option.passwordInputType = "password"; // password | text
setStylesheet(".pasOptionInput {width: 11em;}\n","passwordInputTypeStyle");

merge(config.macros.option.types, {
	'pas': {
		elementType: "input",
		valueField: "value",
		eventName: "onkeyup",
		className: "pasOptionInput",
		typeValue: config.macros.option.passwordInputType,
		create: function(place,type,opt,className,desc) {
			// password field
			config.macros.option.genericCreate(place,'pas',opt,className,desc);
			// checkbox linked with this password "save this password on this computer"
			config.macros.option.genericCreate(place,'chk','chk'+opt,className,desc);			
			// text savePasswordCheckboxLabel
			place.appendChild(document.createTextNode(config.macros.option.passwordCheckboxLabel));
		},
		onChange: config.macros.option.genericOnChange
	}
});

merge(config.optionHandlers['chk'], {
	get: function(name) {
		// is there an option linked with this chk ?
		var opt = name.substr(3);
		if (config.options[opt]) 
			saveOptionCookie(opt);
		return config.options[name] ? "true" : "false";
	}
});

merge(config.optionHandlers, {
	'pas': {
 		get: function(name) {
			if (config.options["chk"+name]) {
				return encodeCookie(config.options[name].toString());
			} else {
				return "";
			}
		},
		set: function(name,value) {config.options[name] = decodeCookie(value);}
	}
});

// need to reload options to load passwordOptions
loadOptionsCookie();

/*
if (!config.options['pasPassword'])
	config.options['pasPassword'] = '';

merge(config.optionsDesc,{
		pasPassword: "Test password"
	});
*/
//}}}
/***
|Name|Plugin: Scientific Notation|
|Created by|BobMcElrath|
|Email|my first name at my last name dot org|
|Location|http://bob.mcelrath.org/tiddlyjsmath-2.0.3.html|
|Version|1.0|
|Requires|[[TiddlyWiki|http://www.tiddlywiki.com]] &ge; 2.0.3, [[jsMath|http://www.math.union.edu/~dpvc/jsMath/]] &ge; 3.0, [[Plugin: jsMath]]|
!Description
This plugin will render numbers expressed in scientific notation, such as {{{3.5483e12}}} using the jsMath plugin to display it in an intuitive way such as 3.5483e12.  You may customize the number of significant figures displayed, as well as "normalize" numbers so that {{{47392.387e9}}} is displayed as 47392.387e9.
!Installation
Install the Requirements, above, add this tiddler to your tiddlywiki, and give it the {{{systemConfig}}} tag.
!History
* 1-Feb-06, version 1.0, Initial release
!Code
***/
//{{{
config.formatters.push({
  name: "scientificNotation",
  match: "\\b[0-9]+\\.[0-9]+[eE][+-]?[0-9]+\\b",
  element: "span",
  className: "math",
  normalize: true,                          // set to 'true' to convert numbers to X.XXX \times 10^{y}
  sigfigs: 3,                               // with this many digits in the mantissa
  handler: function(w) {
    var snRegExp = new RegExp("\\b([0-9]+(?:\\.[0-9]+)?)[eE]([-0-9+]+)\\b");
    var mymatch = snRegExp.exec(w.matchText);
    var mantissa = mymatch[1];
    var exponent = parseInt(mymatch[2]);
    // normalize the number.
    if(this.normalize) {
      mantissa = parseFloat(mantissa);
      while(mantissa > 10.0) {
        mantissa = mantissa / 10.0;
        exponent++; 
      }
      while(mantissa < 1.0) {
        mantissa = mantissa * 10.0;
        exponent--;
      }
      var sigfigsleft = this.sigfigs;
      mantissa = parseInt(mantissa) + "." + (Math.round(Math.pow(10,this.sigfigs-1)*mantissa)+"").substr(1,this.sigfigs-1);
    }
    var e = document.createElement(this.element);
    e.className = this.className;
    if(exponent == 0) {
      e.appendChild(document.createTextNode(mantissa));
    } else {
      e.appendChild(document.createTextNode(mantissa + "\\times 10^{" + exponent + "}"));
    }
    w.output.appendChild(e);
  }
});
//}}}
/***
|Name|Plugin: arXiv Links|
|Created by|BobMcElrath|
|Email|my first name at my last name dot org|
|Location|http://bob.mcelrath.org/tiddlyjsmath-2.0.3.html|
|Version|1.0|
|Requires|[[TiddlyWiki|http://www.tiddlywiki.com]] &ge; 2.0.3|
!Description
This formatting plugin will render links to the [[arXiv|http://www.arxiv.org]] preprint system.  If you type a paper reference such as hep-ph/0509024, it will be rendered as an external link to the abstract of that paper.
!Installation
Add this tiddler to your tiddlywiki, and give it the {{{systemConfig}}} tag.
!History
* 1-Feb-06, version 1.0, Initial release
!Code
***/
//{{{
config.formatters.push({
  name: "arXivLinks",
  match: "\\b(?:astro-ph|cond-mat|hep-ph|hep-th|hep-lat|gr-qc|nucl-ex|nucl-th|quant-ph|(?:cs|math|nlin|physics|q-bio)(?:\\.[A-Z]{2})?)/[0-9]{7}\\b",
  element: "a",
  handler: function(w) {
    var e = createExternalLink(w.output, "http://arxiv.org/abs/"+w.matchText);
    e.target = "_blank"; // open in new window
    w.outputText(e,w.matchStart,w.nextMatch);
  }
});
//}}}
/***
|Name|Plugin: jsMath|
|Created by|BobMcElrath|
|Email|my first name at my last name dot org|
|Location|http://bob.mcelrath.org/tiddlyjsmath.html|
|Version|1.6|
|Requires|[[TiddlyWiki|http://www.tiddlywiki.com]] &ge; 2.0.3, [[jsMath|http://www.math.union.edu/~dpvc/jsMath/]] &ge; 3.0|
!Description
LaTeX is the world standard for specifying, typesetting, and communicating mathematics among scientists, engineers, and mathematicians.  For more information about LaTeX itself, visit the [[LaTeX Project|http://www.latex-project.org/]].  This plugin typesets math using [[jsMath|http://www.math.union.edu/~dpvc/jsMath/]], which is an implementation of the TeX math rules and typesetting in javascript, for your browser.  Notice the small button in the lower right corner which opens its control panel.
!Installation
In addition to this plugin, you must also [[install jsMath|http://www.math.union.edu/~dpvc/jsMath/download/jsMath.html]] on the same server as your TiddlyWiki html file.  If you're using TiddlyWiki without a web server, then the jsMath directory must be placed in the same location as the TiddlyWiki html file.

I also recommend modifying your StyleSheet use serif fonts that are slightly larger than normal, so that the math matches surrounding text, and \\small fonts are not unreadable (as in exponents and subscripts).
{{{
.viewer {
  line-height: 125%;
  font-family: serif;
  font-size: 12pt;
}
}}}

You may also optionally add the following code to load jsMath in [[MarkupPostHead]] if you desire.
{{{
<!--{{{-->
<script src="jsMath/jsMath.js"></script>
<!--}}}-->
}}}
[[Plugin: jsMath]] will normally load jsMath dynamically using AJAX, but adding the above in [[MarkupPostHead]] may be useful if you have jsMath stored in a non-standard location, or if your browser's cross-site origin policy forbids loading files from file URL's using AJAX.  (e.g. Google Chrome)  
!History
* 11-Nov-05, version 1.0, Initial release
* 22-Jan-06, version 1.1, updated for ~TW2.0, tested with jsMath 3.1, editing tiddlywiki.html by hand is no longer necessary.
* 24-Jan-06, version 1.2, fixes for Safari, Konqueror
* 27-Jan-06, version 1.3, improved error handling, detect if ajax was already defined (used by ZiddlyWiki)
* 12-Jul-06, version 1.4, fixed problem with not finding image fonts
* 26-Feb-07, version 1.5, fixed problem with Mozilla "unterminated character class".
* 27-Feb-07, version 1.5.1, Runs compatibly with TW 2.1.0+, by Bram Chen
* 5-May-11, version 1.6, Use a script tag to load in Chrome, use jQuery for ajax
!Examples
|!Source|!Output|h
|{{{The variable $x$ is real.}}}|The variable $x$ is real.|
|{{{The variable \(y\) is complex.}}}|The variable \(y\) is complex.|
|{{{This \[\int_a^b x = \frac{1}{2}(b^2-a^2)\] is an easy integral.}}}|This \[\int_a^b x = \frac{1}{2}(b^2-a^2)\] is an easy integral.|
|{{{This $$\int_a^b \sin x = -(\cos b - \cos a)$$ is another easy integral.}}}|This $$\int_a^b \sin x = -(\cos b - \cos a)$$ is another easy integral.|
|{{{Block formatted equations may also use the 'equation' environment \begin{equation}  \int \tan x = -\ln \cos x \end{equation} }}}|Block formatted equations may also use the 'equation' environment \begin{equation}  \int \tan x = -\ln \cos x \end{equation}|
|{{{Equation arrays are also supported \begin{eqnarray} a &=& b \\ c &=& d \end{eqnarray} }}}|Equation arrays are also supported \begin{eqnarray} a &=& b \\ c &=& d \end{eqnarray} |
|{{{I spent \$7.38 on lunch.}}}|I spent \$7.38 on lunch.|
|{{{I had to insert a backslash (\\) into my document}}}|I had to insert a backslash (\\) into my document|
!Code
***/
//{{{

// Load jsMath
if(typeof jsMath == 'undefined') {
  jsMath = {
    Setup: {inited: 1},          // don't run jsMath.Setup.Body() yet
    Autoload: {root: new String(document.location).replace(/[^\/]*$/,'jsMath/')}  // URL to jsMath directory, change if necessary
  };
  try {
    jQuery.ajax({url: jsMath.Autoload.root+"jsMath.js", dataType: 'script',
                 async: false, error: function(j, s, e) { throw(e); } });
    jsMath.Setup.inited=0;  //  allow jsMath.Setup.Body() to run again
  } catch(e) {
    if(navigator.userAgent.toLowerCase().indexOf('chrome') > -1
        && document.location.toString().indexOf('file') == 0) {
      var jsMathScript = '<!--{{{-->\n<script src="jsMath/jsMath.js"></script>\n<!--}}}-->';
      var MarkupPostHead = store.fetchTiddler("MarkupPostHead");
      if(MarkupPostHead && !MarkupPostHead.match(jsMathScript)) {
        MarkupPostHead.text += jsMathScript;
      } else {
        MarkupPostHead = store.createTiddler("MarkupPostHead");
        MarkupPostHead.text = jsMathScript;
      }
      throw("jsMath added to MarkupPostHead: now save and reload (Google Chrome cross origin file:/// URL workaround)");
    } else {
      alert("jsMath was not found: you must place the 'jsMath' directory in the same place as this file.  "
           +"The error was:\n"+e.name+": "+e.message);
      throw(e);  // abort eval
    }
  }
}

// Define wikifers for latex
config.formatterHelpers.mathFormatHelper = function(w) {
    var e = document.createElement(this.element);
    e.className = this.className;
    var endRegExp = new RegExp(this.terminator, "mg");
    endRegExp.lastIndex = w.matchStart+w.matchLength;
    var matched = endRegExp.exec(w.source);
    if(matched) {
        var txt = w.source.substr(w.matchStart+w.matchLength, 
            matched.index-w.matchStart-w.matchLength);
        if(this.keepdelim) {
          txt = w.source.substr(w.matchStart, matched.index+matched[0].length-w.matchStart);
        }
        e.appendChild(document.createTextNode(txt));
        w.output.appendChild(e);
        w.nextMatch = endRegExp.lastIndex;
    }
}

config.formatters.push({
  name: "displayMath1",
  match: "\\\$\\\$",
  terminator: "\\\$\\\$\\n?", // 2.0 compatability
  termRegExp: "\\\$\\\$\\n?",
  element: "div",
  className: "math",
  handler: config.formatterHelpers.mathFormatHelper
});

config.formatters.push({
  name: "inlineMath1",
  match: "\\\$", 
  terminator: "\\\$", // 2.0 compatability
  termRegExp: "\\\$",
  element: "span",
  className: "math",
  handler: config.formatterHelpers.mathFormatHelper
});

var backslashformatters = new Array(0);

backslashformatters.push({
  name: "inlineMath2",
  match: "\\\\\\\(",
  terminator: "\\\\\\\)", // 2.0 compatability
  termRegExp: "\\\\\\\)",
  element: "span",
  className: "math",
  handler: config.formatterHelpers.mathFormatHelper
});

backslashformatters.push({
  name: "displayMath2",
  match: "\\\\\\\[",
  terminator: "\\\\\\\]\\n?", // 2.0 compatability
  termRegExp: "\\\\\\\]\\n?",
  element: "div",
  className: "math",
  handler: config.formatterHelpers.mathFormatHelper
});

backslashformatters.push({
  name: "displayMath3",
  match: "\\\\begin\\{equation\\}",
  terminator: "\\\\end\\{equation\\}\\n?", // 2.0 compatability
  termRegExp: "\\\\end\\{equation\\}\\n?",
  element: "div",
  className: "math",
  handler: config.formatterHelpers.mathFormatHelper
});

// These can be nested.  e.g. \begin{equation} \begin{array}{ccc} \begin{array}{ccc} ...
backslashformatters.push({
  name: "displayMath4",
  match: "\\\\begin\\{eqnarray\\}",
  terminator: "\\\\end\\{eqnarray\\}\\n?", // 2.0 compatability
  termRegExp: "\\\\end\\{eqnarray\\}\\n?",
  element: "div",
  className: "math",
  keepdelim: true,
  handler: config.formatterHelpers.mathFormatHelper
});

// The escape must come between backslash formatters and regular ones.
// So any latex-like \commands must be added to the beginning of
// backslashformatters here.
backslashformatters.push({
    name: "escape",
    match: "\\\\.",
    handler: function(w) {
        w.output.appendChild(document.createTextNode(w.source.substr(w.matchStart+1,1)));
        w.nextMatch = w.matchStart+2;
    }
});

config.formatters=backslashformatters.concat(config.formatters);

window.wikify = function(source,output,highlightRegExp,tiddler)
{
    if(source && source != "") {
        if(version.major == 2 && version.minor > 0) {
            var wikifier = new Wikifier(source,getParser(tiddler),highlightRegExp,tiddler);
            wikifier.subWikifyUnterm(output);
        } else {
            var wikifier = new Wikifier(source,formatter,highlightRegExp,tiddler);
            wikifier.subWikify(output,null);
        }
        jsMath.ProcessBeforeShowing();
    }
}
//}}}
Consider the following dataset: [img(50%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_original.jpg]]

The correlation matrix and related eigenvalues are given by, [img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_algebra.jpg]]

A one dimensional projection can be obtained by the PCA method, yielding [img(50%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_1d.jpg]]

Including the eigenvector with the two largest eigenvalues, one would have a two dimensional projection:

[img(50%+,aut[img(40%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_2d.jpg]] [img(50%+,aut[img(60%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_2d_flat.jpg]]  

# Develop the high level structure of the algorithm first, and split it in self contained modules.
# Write and test the modules and subroutines/functions individually. Unless the code is very short, do not try to write the entire program at once.
# Use versioning software and repositories for complex code - else keep safety copies. The University operates its own ~GitHub site https://github.umn.edu/ (use your x500 credentials to log in). Keeping versions allows you to go back to previous versions and restore working copies of your project. You can find help at https://help.github.com/articles/set-up-git and https://help.github.com/articles/create-a-repo.
# Stick as much as possible to standard libraries and versions. Many codes need to be run on several different systems and compatibility is a frequent issue.
# Debugging a particular code is the lengthiest part of the development. A code well writen (modular, consistent in notation) makes the task of finding errors simpler. If modules have been tested independently, even better.
# In debugging, if possible, see if the code (or parts) reproduce known solutions.
# Keep in mind that it is possible (although rare) that the computer may make a mistake (typically though internal code optimization).
# Style:
## Indentation is a great way to be able to read, understand and debug the code.
## Comments are useful to document and return to the code at a later time.
## Introduce consisten naming of variables, indices (i,j, ...), loops and arrays.
!!!Data types
|!Type |!FORTRAN descriptor |!C descriptor |!Data size |
|Character | character  | char  | 1 byte |
|Integer | integer  | int | 4 bytes |
|Real (single precision) | real(4) | float | 4 bytes |
|Real (double precision) | real(8) | double | 8 bytes |

Remember: 
* 1 byte = 8 bits
* Modern ~CPUs assume 64 bit words by default. Processor word length is in fact higher. An Intel Haswell processor, for example, has registers that are 512 bit wide (i.e., eight double precision real numbers can be operated on in parallel). This leads to effectively 16 double precision operations per CPU cycle. If interested in details, check http://www.anandtech.com/show/6355/intels-haswell-architecture/8
~PHYS4041. Computational Methods in the Physical Sciences. Fall Semester 2018
University of Minnesota &nbsp; School of Physics and Astronomy <br>
@media print {#mainMenu {display: none ! important;}}
@media print {#topMenu {display: none ! important;}}
@media print {#sidebar {display: none ! important;}}
@media print {#messageArea {display: none ! important;}} 
@media print {#toolbar {display: none ! important;}}
@media print {.header {display: none ! important;}}
@media print {.tiddler .subtitle {display: none ! important;}}
@media print {.tiddler .toolbar {display; none ! important; }}
@media print {.tiddler .tagging {display; none ! important; }}
@media print {.tiddler .tagged {display; none ! important; }}

.justifyright {
  text-align: right;
}
!Course Information
''Instructor:'' Jorge Vi&ntilde;als
''Office:'' Tate Hall, 130-17.
''Contact Information:'' vin&#97;ls&#64;umn.edu
''Office Hours: '' MW, 2:15 p.m. - 3:30 p.m.

''Lecture schedule: '' MTW 1:25 p.m. - 2:15 p.m, Tate Hall B85.
''Lab schedule:'' Th 1:25 p.m. - 2:15 p.m. Keller Hall 1-200.
''TA:'' Brian O'Neill. [[TA course site|http://homepages.spa.umn.edu/~oneill/CompPhys.html]]
''TA Office Hours:'' Tu 2:15 p.m. - 3:30 p.m., Tate Hall 201-02.

!Exams and Final
There will be a take home mid term exam, and one final exam (Monday, December 17, 2018, 1:30 p.m. - 4:30 p.m.). There will be several homework assignments during the semester. Typically, the homework will be due on a Tuesday, one week after being assigned, will involve programming, and may include reading of relevant material. There will be reduced credit for late homework: 75% for up to two days late, 50% for up to on week late. No credit after that. 

!Required Materials

There is no required textbook for this class. Class notes and material posted online should be suffcient. There are, however, a few textbooks that are available online, at the University Library, and at the [[University Bookstore|https://www.bookstores.umn.edu/course/phys-4041-lec-001-computational-methods-twin-cities-fall-2018]] that follow and complement the material in this class.
* F.J. Vesely, "Computational Physics - An Introduction". Kluwer Academic/Plenum Press, New York 2001, second edition. ISBN: 0-306-46631-7. This book is available [[online|https://link-springer-com.ezp2.lib.umn.edu/book/10.1007%2F978-1-4757-2307-6]] for University students.
* R.H. Landau, M. Paez, and C. Bordeianu, "A Survey of Computational Physics". Princeton University Press, New Jersey, 2008. ISBN: 978-0-691-13137-5. On reserve at Walter Library (~QC52.L36 2008).
* B.A. Stickler and E. Schachinger. "Basic Concepts in Computational Physics". Springer, New York (2014), second edition. This book is available [[online|http://link.springer.com/book/10.1007%2F978-3-319-02435-6]] for University students.
A generally good condensed reference, with example codes is:
*  W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery,  [[Numerical recipes|http://numerical.recipes/]]: The Art of Scientific Computing. Cambridge University Press, (2007) third edition. [[Text|http://www.aip.de/groups/soe/local/numres/]].

!Course goals and objectives
This class is an introduction to Scientific Computing from the point of view of Physics applications. It is not a programming class, although some programming experience is assumed and will be necessary during the class. This class is not about learning to use algorithms written by others either, although some reuse, especially of general purpose libraries, will be essential. Rather, the class will focus on the understanding and design of a computing strategy to solve physics problems, the basic numerical analysis that is necessary to develop classes of algorithms, and the implementation of simple but custom application codes.
!!!Topics
# Introduction to scientific computing. Computer architecture, number representation and approximations. Errors and propagation.
# Basic elements of numerical analysis: iteration, local approximation, convergence and accuracy, and stability.
# Brief survey of classical problems: linear and nonlinear systems, interpolation and extrapolation. Differentiation, integration and solution of systems or ordinary differential equations. Boundary and eigenvalue problems.
# Partial differential equations. Methods for elliptic, parabolic and hyperbolic equations. Relaxation methods for boundary value problems.
# Statistical description of data. Probability distributions and their estimation. Data modeling and regression.
# Monte Carlo methods and their quantum extension.
# Elements of data mining. Classification and decision trees. Neural networks. Association analysis in Big Data. Cluster analysis. 
# A high level view of more modern architectures: many core chips, and tightly coupled clusters. Multi threaded and distributed memory programming.

!Grading
| Mid term exam (take home) | 20 % |
| Final exam | 30% |
| Homework assignments | 50 % |

!Attendance policy
Class attendance is required.

!Departamental policies
''ATHLETES'' must provide their official University of Minnesota athletic letter containing the approved competition schedule to their instructor and the staff in office 148. Away exams will be arranged with the athletic adviser traveling with the team. Accommodations will be made for official university sports only (i.e. no accommodations will be made for intramurals, club sports, etc.)
''DISABILITY SERVICES:'' If you have accommodations for this course, please provide the staff in the office (Tate 130) with a copy of your accommodation letter for the current semester. Exams will be arranged according to accommodations and sent to the testing center for administration.

!Mandatory policy information
!!!Student Conduct Code
The University seeks an environment that promotes academic achievement and integrity, that is protective of free inquiry, and that serves the educational mission of the University. Similarly, the University seeks a community that is free from violence, threats, and intimidation; that is respectful of the rights, opportunities, and welfare of students, faculty, staff, and guests of the University; and that does not threaten the physical or mental health or safety of members of the University community. As a student at the University you are expected adhere to Board of Regents Policy: Student Conduct Code. To review the Student Conduct Code, please see: [[Regents Student Conduct Code|http://regents.umn.edu/sites/default/files/policies/Student_Conduct_Code.pdf]]

Note that the conduct code specifically addresses disruptive classroom conduct, which means "engaging in behavior that substantially or repeatedly interrupts either the instructor's ability to teach or student learning. The classroom extends to any setting where a student is engaged in work toward academic credit or satisfaction of program-based requirements or related activities."

!!!Scholastic dishonesty
You are expected to do your own academic work and cite sources as necessary. Failing to do so is scholastic dishonesty. Scholastic dishonesty means plagiarizing; cheating on assignments or examinations; engaging in unauthorized collaboration on academic work; taking, acquiring, or using test materials without faculty permission; submitting false or incomplete records of academic achievement; acting alone or in cooperation with another to falsify records or to obtain dishonestly grades, honors, awards, or professional endorsement; altering, forging, or misusing a University academic record; or fabricating or falsifying data, research procedures, or data analysis ([[Student Conduct Code|http://regents.umn.edu/sites/default/files/policies/Student_Conduct_Code.pdf]]).

If it is determined that a student has cheated, he or she may be given an "F" or an "N" for the course, and may face additional sanctions from the University. For additional information, please see: http://policy.umn.edu/Policies/Education/Education/INSTRUCTORRESP.html.
The Office for Student Conduct and Academic Integrity has compiled a useful list of Frequently Asked Questions pertaining to scholastic dishonesty: http://www1.umn.edu/oscai/integrity/student/index.html. If you have additional questions, please clarify with your instructor for the course. Your instructor can respond to your specific questions regarding what would constitute scholastic dishonesty in the context of a particular class-e.g., whether collaboration on assignments is permitted, requirements and methods for citing sources, if electronic aids are permitted or prohibited during an exam.

!!!Disability Accommodations
The University is committed to providing quality education to all students regardless of ability. Determining appropriate disability accommodations is a collaborative process. You as a student must register with Disability Services and provide documentation of your disability. The course instructor must provide information regarding a course's content, methods, and essential components. The combination of this information will be used by Disability Services to determine appropriate accommodations for a particular student in a particular course. For more information, please reference Disability Services: https://diversity.umn.edu/disability/.

If you have a letter from the Disability Resource Center regarding a special accommodation, please make sure that you share it with the instructor, and the Physics front office.

!!!Use of Personal Electronic Devices in the Classroom
Using personal electronic devices in the classroom setting can hinder instruction and learning, not only for the student using the device but also for other students in the class. To this end, the University establishes the right of each faculty member to determine if and how personal electronic devices are allowed to be used in the classroom. For complete information, please reference: http://www.policy.umn.edu/Policies/Education/Education/STUDENTRESP.html

!!!Makeup Work for Legitimate Absences
Students will not be penalized for absence during the semester due to unavoidable or legitimate circumstances. Such circumstances include verified illness, participation in intercollegiate athletic events, subpoenas, jury duty, military service, bereavement, and religious observances. Such circumstances do not include voting in local, state, or national elections. For complete information, please see: http://policy.umn.edu/Policies/Education/Education/MAKEUPWORK.html.

!!!Appropriate Student Use of Class Notes and Course Materials
Taking notes is a means of recording information but more importantly of personally absorbing and integrating the educational experience. However, broadly disseminating class notes beyond the classroom community or accepting compensation for taking and distributing classroom notes undermines instructor interests in their intellectual work product while not substantially furthering instructor and student interests in effective learning. Such actions violate shared norms and standards of the academic community. For additional information, please see: http://policy.umn.edu/education/studentresp .

!!!Grading and Transcripts
The University utilizes plus and minus grading on a 4.000 cumulative grade point scale in accordance with the following:
|A |4.000 - Represents achievement that is outstanding relative to the level necessary to meet course requirements|
|A- |3.667 |
|B+ |3.333 |
|B |3.000 - Represents achievement that is significantly above the level necessary to meet course requirements |
|B- |2.667 |
|C+ |2.333 |
|C |2.000 - Represents achievement that meets the course requirements in every respect |
|C- |1.667 |
|D+ |1.333 |
|D |1.000 - Represents achievement that is worthy of credit even though it fails to meet fully the course requirements |
|S |Represents achievement that is satisfactory, which is equivalent to a C- or better. |

For additional information, please refer to: http://policy.umn.edu/Policies/Education/Education/GRADINGTRANSCRIPTS.html.

!!!Sexual Harassment
"Sexual harassment" means unwelcome sexual advances, requests for sexual favors, and/or other verbal or physical conduct of a sexual nature. Such conduct has the purpose or effect of unreasonably interfering with an individual's work or academic performance or creating an intimidating, hostile, or offensive working or academic environment in any University activity or program. Such behavior is not acceptable in the University setting. For additional information, please consult Board of Regents Policy: https://policy.umn.edu/hr/sexharassassault .

!!!Equity, Diversity, Equal Opportunity, and Affirmative Action
The University will provide equal access to and opportunity in its programs and facilities, without regard to race, color, creed, religion, national origin, gender, age, marital status, disability, public assistance status, veteran status, sexual orientation, gender identity, or gender expression. For more information, please consult Board of Regents Policy: http://regents.umn.edu/sites/regents.umn.edu/files/policies/Equity_Diversity_EO_AA.pdf .

!!!Mental Health and Stress Management
As a student you may experience a range of issues that can cause barriers to learning, such as strained relationships, increased anxiety, alcohol/drug problems, feeling down, difficulty concentrating and/or lack of motivation. These mental health concerns or stressful events may lead to diminished academic performance and may reduce your ability to participate in daily activities. University of Minnesota services are available to assist you. You can learn more about the broad range of confidential mental health services available on campus via the Student Mental Health Website: http://www.mentalhealth.umn.edu.

!!! Academic Freedom and Responsibility: for courses that do not involve students in research
Academic freedom is a cornerstone of the University. Within the scope and content of the course as defined by the instructor, it includes the freedom to discuss relevant matters in the classroom. Along with this freedom comes responsibility. Students are encouraged to develop the capacity for critical judgment and to engage in a sustained and independent search for truth. Students are free to take reasoned exception to the views offered in any course of study and to reserve judgment about matters of opinion, but they are responsible for learning the content of any course of study for which they are enrolled.

Reports of concerns about academic freedom are taken seriously, and there are individuals and offices available for help. Contact the instructor, the Department Chair, your adviser, the associate dean of the college, or the Vice Provost for Faculty and Academic Affairs in the Office of the Provost.
<<timeline modified 20>>
 <<cloud systemConfig systemTiddlers DiscoveryPackage IconPackage NavigationPackage>>
/***
|Name|TagCloudPlugin|
|Source|http://www.TiddlyTools.com/#TagCloudPlugin|
|Version|1.7.0|
|Author|Eric Shulman|
|Original Author|Clint Checketts|
|License|http://www.TiddlyTools.com/#LegalStatements|
|~CoreVersion|2.1|
|Type|plugin|
|Description|present a 'cloud' of tags (or links) using proportional font display|
!Usage
<<<
{{{
<<cloud type action:... limit:... tag tag tag ...>>
<<cloud type action:... limit:... +TiddlerName>>
<<cloud type action:... limit:... -TiddlerName>>
<<cloud type action:... limit:... =tagvalue>>
}}}
where:
* //type// is a keyword, one of:
** ''tags'' (default) - displays a cloud of tags, based on frequency of use
** ''links'' - displays a cloud of tiddlers, based on number of links //from// each tiddler
** ''references'' - displays a cloud of tiddlers, based on number of links //to// each tiddler
* ''action:popup'' (default) - clicking a cloud item shows a popup with links to related tiddlers<br>//or//<br> ''action:goto'' - clicking a cloud item immediately opens the tiddler corresponding to that item
* ''limit:N'' (optional) - restricts the cloud display to only show the N most popular tags/links
* ''tag tag tag...'' (or ''title title title'' if ''links''/''references'' is used)<br>shows all tags/links in the document //except// for those listed as macro parameters
* ''+TiddlerName''<br>show only tags/links read from a space-separated, bracketed list stored in a separate tiddler.
* ''-TiddlerName''<br>show all tags/links //except// those read from a space-separated, bracketed list stored in a separate tiddler.
* ''=tagvalue'' (//only if type=''tags''//)<br>shows only tags that are themselves tagged with the indicated tag value (i.e., ~TagglyTagging usage)
//note: for backward-compatibility, you can also use the macro {{{<<tagCloud ...>>}}} in place of {{{<<cloud ...>>}}}//
<<<
!Examples
<<<
//all tags excluding<<tag systemConfig>>, <<tag excludeMissing>> and <<tag script>>//
{{{<<cloud systemConfig excludeMissing script>>}}}
{{groupbox{<<cloud systemConfig excludeMissing script>>}}}
//top 10 tags excluding<<tag systemConfig>>, <<tag excludeMissing>> and <<tag script>>//
{{{<<cloud limit:10 systemConfig excludeMissing script>>}}}
{{groupbox{<<cloud limit:10 systemConfig excludeMissing script>>}}}
//tags listed in// [[FavoriteTags]]
{{{<<cloud +FavoriteTags>>}}}
{{groupbox{<<cloud +FavoriteTags>>}}}
//tags NOT listed in// [[FavoriteTags]]
{{{<<cloud -FavoriteTags>>}}}
{{groupbox{<<cloud -FavoriteTags>>}}}
//links to tiddlers tagged with 'package'//
{{{<<cloud action:goto =package>>}}}
{{groupbox{<<cloud action:goto =package>>}}}
//top 20 most referenced tiddlers//
{{{<<cloud references limit:20>>}}}
{{groupbox{<<cloud references limit:20>>}}}
//top 20 tiddlers that contain the most links//
{{{<<cloud links limit:20>>}}}
{{groupbox{<<cloud links limit:20>>}}}
<<<
!Revisions
<<<
2009.07.17 [1.7.0] added {{{-TiddlerName}}} parameter to exclude tags that are listed in the indicated tiddler
2009.02.26 [1.6.0] added {{{action:...}}} parameter to apply popup vs. goto action when clicking cloud items
2009.02.05 [1.5.0] added ability to show links or back-links (references) instead of tags and renamed macro to {{{<<cloud>>}}} to reflect more generalized usage.
2008.12.16 [1.4.2] corrected group calculation to prevent 'group=0' error
2008.12.16 [1.4.1] revised tag filtering so excluded tags don't affect calculations
2008.12.15 [1.4.0] added {{{limit:...}}} parameter to restrict the number of tags displayed to the top N most popular
2008.11.15 [1.3.0] added {{{+TiddlerName}}} parameter to include only tags that are listed in the indicated tiddler
2008.09.05 [1.2.0] added '=tagname' parameter to include only tags that are themselves tagged with the specified value (i.e., ~TagglyTagging usage)
2008.07.03 [1.1.0] added 'segments' property to macro object.  Extensive code cleanup
<<<
!Code
***/
//{{{
version.extensions.TagCloudPlugin= {major: 1, minor: 7 , revision: 0, date: new Date(2009,7,17)};
//Originally created by Clint Checketts, contributions by Jonny Leroy and Eric Shulman
//Currently maintained and enhanced by Eric Shulman
//}}}
//{{{
config.macros.cloud = {
	tagstip: "%1 tiddlers tagged with '%0'",
	refslabel: " (%0 references)",
	refstip: "%1 tiddlers have links to '%0'",
	linkslabel: " (%0 links)",
	linkstip: "'%0' has links to %1 other tiddlers",
	groups: 9,
	init: function() {
		config.macros.tagCloud=config.macros.cloud; // for backward-compatibility
		config.shadowTiddlers.TagCloud='<<cloud>>';
		config.shadowTiddlers.StyleSheetTagCloud=
			'/*{{{*/\n'
			+'.tagCloud span {line-height: 3.5em; margin:3px;}\n'
			+'.tagCloud1{font-size: 80%;}\n'
			+'.tagCloud2{font-size: 100%;}\n'
			+'.tagCloud3{font-size: 120%;}\n'
			+'.tagCloud4{font-size: 140%;}\n'
			+'.tagCloud5{font-size: 160%;}\n'
			+'.tagCloud6{font-size: 180%;}\n'
			+'.tagCloud7{font-size: 200%;}\n'
			+'.tagCloud8{font-size: 220%;}\n'
			+'.tagCloud9{font-size: 240%;}\n'
			+'/*}}}*/\n';
		setStylesheet(store.getTiddlerText('StyleSheetTagCloud'),'tagCloudsStyles');
	},
	getLinks: function(tiddler) { // get list of links to existing tiddlers and shadows
		if (!tiddler.linksUpdated) tiddler.changed();
		var list=[]; for (var i=0; i<tiddler.links.length; i++) {
			var title=tiddler.links[i];
			if (store.isShadowTiddler(title)||store.tiddlerExists(title))
				list.push(title);
		}
		return list;
	},
	handler: function(place,macroName,params) {
		// unpack params
		var inc=[]; var ex=[]; var limit=0; var action='popup';
		var links=(params[0]&&params[0].toLowerCase()=='links'); if (links) params.shift();
		var refs=(params[0]&&params[0].toLowerCase()=='references'); if (refs) params.shift();
		if (params[0]&&params[0].substr(0,7).toLowerCase()=='action:')
			action=params.shift().substr(7).toLowerCase();
		if (params[0]&&params[0].substr(0,6).toLowerCase()=='limit:')
			limit=parseInt(params.shift().substr(6));
		while (params.length) {
			if (params[0].substr(0,1)=='+') { // read taglist from tiddler
				inc=inc.concat(store.getTiddlerText(params[0].substr(1),'').readBracketedList());
			} else if (params[0].substr(0,1)=='-') { // exclude taglist from tiddler
				ex=ex.concat(store.getTiddlerText(params[0].substr(1),'').readBracketedList());
			} else if (params[0].substr(0,1)=='=') { // get tag list using tagged tags
				var tagged=store.getTaggedTiddlers(params[0].substr(1));
				for (var t=0; t<tagged.length; t++) inc.push(tagged[t].title);
			} else ex.push(params[0]); // exclude params
			params.shift();
		}
		// get all items, include/exclude specific items
		var items=[];
		var list=(links||refs)?store.getTiddlers('title','excludeLists'):store.getTags();
		for (var t=0; t<list.length; t++) {
			var title=(links||refs)?list[t].title:list[t][0];
			if (links)	var count=this.getLinks(list[t]).length;
			else if (refs)	var count=store.getReferringTiddlers(title).length;
			else 		var count=list[t][1];
			if ((!inc.length||inc.contains(title))&&(!ex.length||!ex.contains(title)))
				items.push({ title:title, count:count });
		}
		if(!items.length) return;
		// sort by decending count, limit results (optional)
		items=items.sort(function(a,b){return(a.count==b.count)?0:(a.count>b.count?-1:1);});
		while (limit && items.length>limit) items.pop();
		// find min/max and group size
		var most=items[0].count;
		var least=items[items.length-1].count;
		var groupSize=(most-least+1)/this.groups;
		// sort by title and draw the cloud of items
		items=items.sort(function(a,b){return(a.title==b.title)?0:(a.title>b.title?1:-1);});
		var cloudWrapper = createTiddlyElement(place,'div',null,'tagCloud',null);
		for (var t=0; t<items.length; t++) {
			cloudWrapper.appendChild(document.createTextNode(' '));
			var group=Math.ceil((items[t].count-least)/groupSize)||1;
			var className='tagCloudtag tagCloud'+group;
			var tip=refs?this.refstip:links?this.linkstip:this.tagstip;
			tip=tip.format([items[t].title,items[t].count]);
			if (action=='goto') { // TAG/LINK/REFERENCES GOTO
				var btn=createTiddlyLink(cloudWrapper,items[t].title,true,className);
				btn.title=tip;
				btn.style.fontWeight='normal';
			} else if (!links&&!refs) { // TAG POPUP
				var btn=createTiddlyButton(cloudWrapper,items[t].title,tip,onClickTag,className);
				btn.setAttribute('tag',items[t].title);
			} else { // LINK/REFERENCES POPUP
				var btn=createTiddlyButton(cloudWrapper,items[t].title,tip,
					function(ev) { var e=ev||window.event; var cmt=config.macros.cloud;
						var popup = Popup.create(this);
						var title = this.getAttribute('tiddler');
						var count = this.getAttribute('count');
						var refs  = this.getAttribute('refs')=='T';
						var links = this.getAttribute('links')=='T';
						var label = (refs?cmt.refslabel:cmt.linkslabel).format([count]);
						createTiddlyLink(popup,title,true);
						createTiddlyText(popup,label);
						createTiddlyElement(popup,'hr');
						if (refs) {
							popup.setAttribute('tiddler',title);
							config.commands.references.handlePopup(popup,title);
						}
						if (links) {
							var tiddler = store.fetchTiddler(title);
							var links=config.macros.cloud.getLinks(tiddler);
							for(var i=0;i<links.length;i++)
								createTiddlyLink(createTiddlyElement(popup,'li'),
									links[i],true);
						}
						Popup.show();
						e.cancelBubble=true; if(e.stopPropagation) e.stopPropagation();
						return false;
					}, className);
				btn.setAttribute('tiddler',items[t].title);
				btn.setAttribute('count',items[t].count);
				btn.setAttribute('refs',refs?'T':'F');
				btn.setAttribute('links',links?'T':'F');
				btn.title=tip;
			}
		}
	}
};
//}}}
/***
Description: Contains the stuff you need to use Tiddlyspot
Note, you also need UploadPlugin, PasswordOptionPlugin and LoadRemoteFileThroughProxy
from http://tiddlywiki.bidix.info for a complete working Tiddlyspot site.
***/
//{{{

// edit this if you are migrating sites or retrofitting an existing TW
config.tiddlyspotSiteId = 'phys4041';

// make it so you can by default see edit controls via http
config.options.chkHttpReadOnly = false;
window.readOnly = false; // make sure of it (for tw 2.2)
window.showBackstage = true; // show backstage too

// disable autosave in d3
if (window.location.protocol != "file:")
	config.options.chkGTDLazyAutoSave = false;

// tweak shadow tiddlers to add upload button, password entry box etc
with (config.shadowTiddlers) {
	SiteUrl = 'http://'+config.tiddlyspotSiteId+'.tiddlyspot.com';
	SideBarOptions = SideBarOptions.replace(/(<<saveChanges>>)/,"$1<<tiddler TspotSidebar>>");
	OptionsPanel = OptionsPanel.replace(/^/,"<<tiddler TspotOptions>>");
	DefaultTiddlers = DefaultTiddlers.replace(/^/,"[[WelcomeToTiddlyspot]] ");
	MainMenu = MainMenu.replace(/^/,"[[WelcomeToTiddlyspot]] ");
}

// create some shadow tiddler content
merge(config.shadowTiddlers,{

'TspotControls':[
 "| tiddlyspot password:|<<option pasUploadPassword>>|",
 "| site management:|<<upload http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/store.cgi index.html . .  " + config.tiddlyspotSiteId + ">>//(requires tiddlyspot password)//<br>[[control panel|http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/controlpanel]], [[download (go offline)|http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/download]]|",
 "| links:|[[tiddlyspot.com|http://tiddlyspot.com/]], [[FAQs|http://faq.tiddlyspot.com/]], [[blog|http://tiddlyspot.blogspot.com/]], email [[support|mailto:support@tiddlyspot.com]] & [[feedback|mailto:feedback@tiddlyspot.com]], [[donate|http://tiddlyspot.com/?page=donate]]|"
].join("\n"),

'TspotOptions':[
 "tiddlyspot password:",
 "<<option pasUploadPassword>>",
 ""
].join("\n"),

'TspotSidebar':[
 "<<upload http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/store.cgi index.html . .  " + config.tiddlyspotSiteId + ">><html><a href='http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/download' class='button'>download</a></html>"
].join("\n"),

'WelcomeToTiddlyspot':[
 "This document is a ~TiddlyWiki from tiddlyspot.com.  A ~TiddlyWiki is an electronic notebook that is great for managing todo lists, personal information, and all sorts of things.",
 "",
 "@@font-weight:bold;font-size:1.3em;color:#444; //What now?// &nbsp;&nbsp;@@ Before you can save any changes, you need to enter your password in the form below.  Then configure privacy and other site settings at your [[control panel|http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/controlpanel]] (your control panel username is //" + config.tiddlyspotSiteId + "//).",
 "<<tiddler TspotControls>>",
 "See also GettingStarted.",
 "",
 "@@font-weight:bold;font-size:1.3em;color:#444; //Working online// &nbsp;&nbsp;@@ You can edit this ~TiddlyWiki right now, and save your changes using the \"save to web\" button in the column on the right.",
 "",
 "@@font-weight:bold;font-size:1.3em;color:#444; //Working offline// &nbsp;&nbsp;@@ A fully functioning copy of this ~TiddlyWiki can be saved onto your hard drive or USB stick.  You can make changes and save them locally without being connected to the Internet.  When you're ready to sync up again, just click \"upload\" and your ~TiddlyWiki will be saved back to tiddlyspot.com.",
 "",
 "@@font-weight:bold;font-size:1.3em;color:#444; //Help!// &nbsp;&nbsp;@@ Find out more about ~TiddlyWiki at [[TiddlyWiki.com|http://tiddlywiki.com]].  Also visit [[TiddlyWiki.org|http://tiddlywiki.org]] for documentation on learning and using ~TiddlyWiki. New users are especially welcome on the [[TiddlyWiki mailing list|http://groups.google.com/group/TiddlyWiki]], which is an excellent place to ask questions and get help.  If you have a tiddlyspot related problem email [[tiddlyspot support|mailto:support@tiddlyspot.com]].",
 "",
 "@@font-weight:bold;font-size:1.3em;color:#444; //Enjoy :)// &nbsp;&nbsp;@@ We hope you like using your tiddlyspot.com site.  Please email [[feedback@tiddlyspot.com|mailto:feedback@tiddlyspot.com]] with any comments or suggestions."
].join("\n")

});
//}}}
| !date | !user | !location | !storeUrl | !uploadDir | !toFilename | !backupdir | !origin |
| 31/10/2018 09:10:47 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
| 02/11/2018 07:49:36 | Brian | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | failed |
| 02/11/2018 07:49:47 | Brian | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | failed |
| 02/11/2018 07:50:00 | Brian | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | failed |
| 02/11/2018 07:50:22 | Brian | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
| 06/11/2018 15:41:06 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
| 07/11/2018 21:10:47 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | ok |
| 07/11/2018 21:11:29 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
| 08/11/2018 13:24:40 | Brian | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
| 11/11/2018 21:24:47 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
/***
|''Name:''|UploadPlugin|
|''Description:''|Save to web a TiddlyWiki|
|''Version:''|4.1.3|
|''Date:''|Feb 24, 2008|
|''Source:''|http://tiddlywiki.bidix.info/#UploadPlugin|
|''Documentation:''|http://tiddlywiki.bidix.info/#UploadPluginDoc|
|''Author:''|BidiX (BidiX (at) bidix (dot) info)|
|''License:''|[[BSD open source license|http://tiddlywiki.bidix.info/#%5B%5BBSD%20open%20source%20license%5D%5D ]]|
|''~CoreVersion:''|2.2.0|
|''Requires:''|PasswordOptionPlugin|
***/
//{{{
version.extensions.UploadPlugin = {
	major: 4, minor: 1, revision: 3,
	date: new Date("Feb 24, 2008"),
	source: 'http://tiddlywiki.bidix.info/#UploadPlugin',
	author: 'BidiX (BidiX (at) bidix (dot) info',
	coreVersion: '2.2.0'
};

//
// Environment
//

if (!window.bidix) window.bidix = {}; // bidix namespace
bidix.debugMode = false;	// true to activate both in Plugin and UploadService
	
//
// Upload Macro
//

config.macros.upload = {
// default values
	defaultBackupDir: '',	//no backup
	defaultStoreScript: "store.php",
	defaultToFilename: "index.html",
	defaultUploadDir: ".",
	authenticateUser: true	// UploadService Authenticate User
};
	
config.macros.upload.label = {
	promptOption: "Save and Upload this TiddlyWiki with UploadOptions",
	promptParamMacro: "Save and Upload this TiddlyWiki in %0",
	saveLabel: "save to web", 
	saveToDisk: "save to disk",
	uploadLabel: "upload"	
};

config.macros.upload.messages = {
	noStoreUrl: "No store URL in parmeters or options",
	usernameOrPasswordMissing: "Username or password missing"
};

config.macros.upload.handler = function(place,macroName,params) {
	if (readOnly)
		return;
	var label;
	if (document.location.toString().substr(0,4) == "http") 
		label = this.label.saveLabel;
	else
		label = this.label.uploadLabel;
	var prompt;
	if (params[0]) {
		prompt = this.label.promptParamMacro.toString().format([this.destFile(params[0], 
			(params[1] ? params[1]:bidix.basename(window.location.toString())), params[3])]);
	} else {
		prompt = this.label.promptOption;
	}
	createTiddlyButton(place, label, prompt, function() {config.macros.upload.action(params);}, null, null, this.accessKey);
};

config.macros.upload.action = function(params)
{
		// for missing macro parameter set value from options
		if (!params) params = {};
		var storeUrl = params[0] ? params[0] : config.options.txtUploadStoreUrl;
		var toFilename = params[1] ? params[1] : config.options.txtUploadFilename;
		var backupDir = params[2] ? params[2] : config.options.txtUploadBackupDir;
		var uploadDir = params[3] ? params[3] : config.options.txtUploadDir;
		var username = params[4] ? params[4] : config.options.txtUploadUserName;
		var password = config.options.pasUploadPassword; // for security reason no password as macro parameter	
		// for still missing parameter set default value
		if ((!storeUrl) && (document.location.toString().substr(0,4) == "http")) 
			storeUrl = bidix.dirname(document.location.toString())+'/'+config.macros.upload.defaultStoreScript;
		if (storeUrl.substr(0,4) != "http")
			storeUrl = bidix.dirname(document.location.toString()) +'/'+ storeUrl;
		if (!toFilename)
			toFilename = bidix.basename(window.location.toString());
		if (!toFilename)
			toFilename = config.macros.upload.defaultToFilename;
		if (!uploadDir)
			uploadDir = config.macros.upload.defaultUploadDir;
		if (!backupDir)
			backupDir = config.macros.upload.defaultBackupDir;
		// report error if still missing
		if (!storeUrl) {
			alert(config.macros.upload.messages.noStoreUrl);
			clearMessage();
			return false;
		}
		if (config.macros.upload.authenticateUser && (!username || !password)) {
			alert(config.macros.upload.messages.usernameOrPasswordMissing);
			clearMessage();
			return false;
		}
		bidix.upload.uploadChanges(false,null,storeUrl, toFilename, uploadDir, backupDir, username, password); 
		return false; 
};

config.macros.upload.destFile = function(storeUrl, toFilename, uploadDir) 
{
	if (!storeUrl)
		return null;
		var dest = bidix.dirname(storeUrl);
		if (uploadDir && uploadDir != '.')
			dest = dest + '/' + uploadDir;
		dest = dest + '/' + toFilename;
	return dest;
};

//
// uploadOptions Macro
//

config.macros.uploadOptions = {
	handler: function(place,macroName,params) {
		var wizard = new Wizard();
		wizard.createWizard(place,this.wizardTitle);
		wizard.addStep(this.step1Title,this.step1Html);
		var markList = wizard.getElement("markList");
		var listWrapper = document.createElement("div");
		markList.parentNode.insertBefore(listWrapper,markList);
		wizard.setValue("listWrapper",listWrapper);
		this.refreshOptions(listWrapper,false);
		var uploadCaption;
		if (document.location.toString().substr(0,4) == "http") 
			uploadCaption = config.macros.upload.label.saveLabel;
		else
			uploadCaption = config.macros.upload.label.uploadLabel;
		
		wizard.setButtons([
				{caption: uploadCaption, tooltip: config.macros.upload.label.promptOption, 
					onClick: config.macros.upload.action},
				{caption: this.cancelButton, tooltip: this.cancelButtonPrompt, onClick: this.onCancel}
				
			]);
	},
	options: [
		"txtUploadUserName",
		"pasUploadPassword",
		"txtUploadStoreUrl",
		"txtUploadDir",
		"txtUploadFilename",
		"txtUploadBackupDir",
		"chkUploadLog",
		"txtUploadLogMaxLine"		
	],
	refreshOptions: function(listWrapper) {
		var opts = [];
		for(i=0; i<this.options.length; i++) {
			var opt = {};
			opts.push();
			opt.option = "";
			n = this.options[i];
			opt.name = n;
			opt.lowlight = !config.optionsDesc[n];
			opt.description = opt.lowlight ? this.unknownDescription : config.optionsDesc[n];
			opts.push(opt);
		}
		var listview = ListView.create(listWrapper,opts,this.listViewTemplate);
		for(n=0; n<opts.length; n++) {
			var type = opts[n].name.substr(0,3);
			var h = config.macros.option.types[type];
			if (h && h.create) {
				h.create(opts[n].colElements['option'],type,opts[n].name,opts[n].name,"no");
			}
		}
		
	},
	onCancel: function(e)
	{
		backstage.switchTab(null);
		return false;
	},
	
	wizardTitle: "Upload with options",
	step1Title: "These options are saved in cookies in your browser",
	step1Html: "<input type='hidden' name='markList'></input><br>",
	cancelButton: "Cancel",
	cancelButtonPrompt: "Cancel prompt",
	listViewTemplate: {
		columns: [
			{name: 'Description', field: 'description', title: "Description", type: 'WikiText'},
			{name: 'Option', field: 'option', title: "Option", type: 'String'},
			{name: 'Name', field: 'name', title: "Name", type: 'String'}
			],
		rowClasses: [
			{className: 'lowlight', field: 'lowlight'} 
			]}
};

//
// upload functions
//

if (!bidix.upload) bidix.upload = {};

if (!bidix.upload.messages) bidix.upload.messages = {
	//from saving
	invalidFileError: "The original file '%0' does not appear to be a valid TiddlyWiki",
	backupSaved: "Backup saved",
	backupFailed: "Failed to upload backup file",
	rssSaved: "RSS feed uploaded",
	rssFailed: "Failed to upload RSS feed file",
	emptySaved: "Empty template uploaded",
	emptyFailed: "Failed to upload empty template file",
	mainSaved: "Main TiddlyWiki file uploaded",
	mainFailed: "Failed to upload main TiddlyWiki file. Your changes have not been saved",
	//specific upload
	loadOriginalHttpPostError: "Can't get original file",
	aboutToSaveOnHttpPost: 'About to upload on %0 ...',
	storePhpNotFound: "The store script '%0' was not found."
};

bidix.upload.uploadChanges = function(onlyIfDirty,tiddlers,storeUrl,toFilename,uploadDir,backupDir,username,password)
{
	var callback = function(status,uploadParams,original,url,xhr) {
		if (!status) {
			displayMessage(bidix.upload.messages.loadOriginalHttpPostError);
			return;
		}
		if (bidix.debugMode) 
			alert(original.substr(0,500)+"\n...");
		// Locate the storeArea div's 
		var posDiv = locateStoreArea(original);
		if((posDiv[0] == -1) || (posDiv[1] == -1)) {
			alert(config.messages.invalidFileError.format([localPath]));
			return;
		}
		bidix.upload.uploadRss(uploadParams,original,posDiv);
	};
	
	if(onlyIfDirty && !store.isDirty())
		return;
	clearMessage();
	// save on localdisk ?
	if (document.location.toString().substr(0,4) == "file") {
		var path = document.location.toString();
		var localPath = getLocalPath(path);
		saveChanges();
	}
	// get original
	var uploadParams = new Array(storeUrl,toFilename,uploadDir,backupDir,username,password);
	var originalPath = document.location.toString();
	// If url is a directory : add index.html
	if (originalPath.charAt(originalPath.length-1) == "/")
		originalPath = originalPath + "index.html";
	var dest = config.macros.upload.destFile(storeUrl,toFilename,uploadDir);
	var log = new bidix.UploadLog();
	log.startUpload(storeUrl, dest, uploadDir,  backupDir);
	displayMessage(bidix.upload.messages.aboutToSaveOnHttpPost.format([dest]));
	if (bidix.debugMode) 
		alert("about to execute Http - GET on "+originalPath);
	var r = doHttp("GET",originalPath,null,null,username,password,callback,uploadParams,null);
	if (typeof r == "string")
		displayMessage(r);
	return r;
};

bidix.upload.uploadRss = function(uploadParams,original,posDiv) 
{
	var callback = function(status,params,responseText,url,xhr) {
		if(status) {
			var destfile = responseText.substring(responseText.indexOf("destfile:")+9,responseText.indexOf("\n", responseText.indexOf("destfile:")));
			displayMessage(bidix.upload.messages.rssSaved,bidix.dirname(url)+'/'+destfile);
			bidix.upload.uploadMain(params[0],params[1],params[2]);
		} else {
			displayMessage(bidix.upload.messages.rssFailed);			
		}
	};
	// do uploadRss
	if(config.options.chkGenerateAnRssFeed) {
		var rssPath = uploadParams[1].substr(0,uploadParams[1].lastIndexOf(".")) + ".xml";
		var rssUploadParams = new Array(uploadParams[0],rssPath,uploadParams[2],'',uploadParams[4],uploadParams[5]);
		var rssString = generateRss();
		// no UnicodeToUTF8 conversion needed when location is "file" !!!
		if (document.location.toString().substr(0,4) != "file")
			rssString = convertUnicodeToUTF8(rssString);	
		bidix.upload.httpUpload(rssUploadParams,rssString,callback,Array(uploadParams,original,posDiv));
	} else {
		bidix.upload.uploadMain(uploadParams,original,posDiv);
	}
};

bidix.upload.uploadMain = function(uploadParams,original,posDiv) 
{
	var callback = function(status,params,responseText,url,xhr) {
		var log = new bidix.UploadLog();
		if(status) {
			// if backupDir specified
			if ((params[3]) && (responseText.indexOf("backupfile:") > -1))  {
				var backupfile = responseText.substring(responseText.indexOf("backupfile:")+11,responseText.indexOf("\n", responseText.indexOf("backupfile:")));
				displayMessage(bidix.upload.messages.backupSaved,bidix.dirname(url)+'/'+backupfile);
			}
			var destfile = responseText.substring(responseText.indexOf("destfile:")+9,responseText.indexOf("\n", responseText.indexOf("destfile:")));
			displayMessage(bidix.upload.messages.mainSaved,bidix.dirname(url)+'/'+destfile);
			store.setDirty(false);
			log.endUpload("ok");
		} else {
			alert(bidix.upload.messages.mainFailed);
			displayMessage(bidix.upload.messages.mainFailed);
			log.endUpload("failed");			
		}
	};
	// do uploadMain
	var revised = bidix.upload.updateOriginal(original,posDiv);
	bidix.upload.httpUpload(uploadParams,revised,callback,uploadParams);
};

bidix.upload.httpUpload = function(uploadParams,data,callback,params)
{
	var localCallback = function(status,params,responseText,url,xhr) {
		url = (url.indexOf("nocache=") < 0 ? url : url.substring(0,url.indexOf("nocache=")-1));
		if (xhr.status == 404)
			alert(bidix.upload.messages.storePhpNotFound.format([url]));
		if ((bidix.debugMode) || (responseText.indexOf("Debug mode") >= 0 )) {
			alert(responseText);
			if (responseText.indexOf("Debug mode") >= 0 )
				responseText = responseText.substring(responseText.indexOf("\n\n")+2);
		} else if (responseText.charAt(0) != '0') 
			alert(responseText);
		if (responseText.charAt(0) != '0')
			status = null;
		callback(status,params,responseText,url,xhr);
	};
	// do httpUpload
	var boundary = "---------------------------"+"AaB03x";	
	var uploadFormName = "UploadPlugin";
	// compose headers data
	var sheader = "";
	sheader += "--" + boundary + "\r\nContent-disposition: form-data; name=\"";
	sheader += uploadFormName +"\"\r\n\r\n";
	sheader += "backupDir="+uploadParams[3] +
				";user=" + uploadParams[4] +
				";password=" + uploadParams[5] +
				";uploaddir=" + uploadParams[2];
	if (bidix.debugMode)
		sheader += ";debug=1";
	sheader += ";;\r\n"; 
	sheader += "\r\n" + "--" + boundary + "\r\n";
	sheader += "Content-disposition: form-data; name=\"userfile\"; filename=\""+uploadParams[1]+"\"\r\n";
	sheader += "Content-Type: text/html;charset=UTF-8" + "\r\n";
	sheader += "Content-Length: " + data.length + "\r\n\r\n";
	// compose trailer data
	var strailer = new String();
	strailer = "\r\n--" + boundary + "--\r\n";
	data = sheader + data + strailer;
	if (bidix.debugMode) alert("about to execute Http - POST on "+uploadParams[0]+"\n with \n"+data.substr(0,500)+ " ... ");
	var r = doHttp("POST",uploadParams[0],data,"multipart/form-data; ;charset=UTF-8; boundary="+boundary,uploadParams[4],uploadParams[5],localCallback,params,null);
	if (typeof r == "string")
		displayMessage(r);
	return r;
};

// same as Saving's updateOriginal but without convertUnicodeToUTF8 calls
bidix.upload.updateOriginal = function(original, posDiv)
{
	if (!posDiv)
		posDiv = locateStoreArea(original);
	if((posDiv[0] == -1) || (posDiv[1] == -1)) {
		alert(config.messages.invalidFileError.format([localPath]));
		return;
	}
	var revised = original.substr(0,posDiv[0] + startSaveArea.length) + "\n" +
				store.allTiddlersAsHtml() + "\n" +
				original.substr(posDiv[1]);
	var newSiteTitle = getPageTitle().htmlEncode();
	revised = revised.replaceChunk("<title"+">","</title"+">"," " + newSiteTitle + " ");
	revised = updateMarkupBlock(revised,"PRE-HEAD","MarkupPreHead");
	revised = updateMarkupBlock(revised,"POST-HEAD","MarkupPostHead");
	revised = updateMarkupBlock(revised,"PRE-BODY","MarkupPreBody");
	revised = updateMarkupBlock(revised,"POST-SCRIPT","MarkupPostBody");
	return revised;
};

//
// UploadLog
// 
// config.options.chkUploadLog :
//		false : no logging
//		true : logging
// config.options.txtUploadLogMaxLine :
//		-1 : no limit
//      0 :  no Log lines but UploadLog is still in place
//		n :  the last n lines are only kept
//		NaN : no limit (-1)

bidix.UploadLog = function() {
	if (!config.options.chkUploadLog) 
		return; // this.tiddler = null
	this.tiddler = store.getTiddler("UploadLog");
	if (!this.tiddler) {
		this.tiddler = new Tiddler();
		this.tiddler.title = "UploadLog";
		this.tiddler.text = "| !date | !user | !location | !storeUrl | !uploadDir | !toFilename | !backupdir | !origin |";
		this.tiddler.created = new Date();
		this.tiddler.modifier = config.options.txtUserName;
		this.tiddler.modified = new Date();
		store.addTiddler(this.tiddler);
	}
	return this;
};

bidix.UploadLog.prototype.addText = function(text) {
	if (!this.tiddler)
		return;
	// retrieve maxLine when we need it
	var maxLine = parseInt(config.options.txtUploadLogMaxLine,10);
	if (isNaN(maxLine))
		maxLine = -1;
	// add text
	if (maxLine != 0) 
		this.tiddler.text = this.tiddler.text + text;
	// Trunck to maxLine
	if (maxLine >= 0) {
		var textArray = this.tiddler.text.split('\n');
		if (textArray.length > maxLine + 1)
			textArray.splice(1,textArray.length-1-maxLine);
			this.tiddler.text = textArray.join('\n');		
	}
	// update tiddler fields
	this.tiddler.modifier = config.options.txtUserName;
	this.tiddler.modified = new Date();
	store.addTiddler(this.tiddler);
	// refresh and notifiy for immediate update
	story.refreshTiddler(this.tiddler.title);
	store.notify(this.tiddler.title, true);
};

bidix.UploadLog.prototype.startUpload = function(storeUrl, toFilename, uploadDir,  backupDir) {
	if (!this.tiddler)
		return;
	var now = new Date();
	var text = "\n| ";
	var filename = bidix.basename(document.location.toString());
	if (!filename) filename = '/';
	text += now.formatString("0DD/0MM/YYYY 0hh:0mm:0ss") +" | ";
	text += config.options.txtUserName + " | ";
	text += "[["+filename+"|"+location + "]] |";
	text += " [[" + bidix.basename(storeUrl) + "|" + storeUrl + "]] | ";
	text += uploadDir + " | ";
	text += "[[" + bidix.basename(toFilename) + " | " +toFilename + "]] | ";
	text += backupDir + " |";
	this.addText(text);
};

bidix.UploadLog.prototype.endUpload = function(status) {
	if (!this.tiddler)
		return;
	this.addText(" "+status+" |");
};

//
// Utilities
// 

bidix.checkPlugin = function(plugin, major, minor, revision) {
	var ext = version.extensions[plugin];
	if (!
		(ext  && 
			((ext.major > major) || 
			((ext.major == major) && (ext.minor > minor))  ||
			((ext.major == major) && (ext.minor == minor) && (ext.revision >= revision))))) {
			// write error in PluginManager
			if (pluginInfo)
				pluginInfo.log.push("Requires " + plugin + " " + major + "." + minor + "." + revision);
			eval(plugin); // generate an error : "Error: ReferenceError: xxxx is not defined"
	}
};

bidix.dirname = function(filePath) {
	if (!filePath) 
		return;
	var lastpos;
	if ((lastpos = filePath.lastIndexOf("/")) != -1) {
		return filePath.substring(0, lastpos);
	} else {
		return filePath.substring(0, filePath.lastIndexOf("\\"));
	}
};

bidix.basename = function(filePath) {
	if (!filePath) 
		return;
	var lastpos;
	if ((lastpos = filePath.lastIndexOf("#")) != -1) 
		filePath = filePath.substring(0, lastpos);
	if ((lastpos = filePath.lastIndexOf("/")) != -1) {
		return filePath.substring(lastpos + 1);
	} else
		return filePath.substring(filePath.lastIndexOf("\\")+1);
};

bidix.initOption = function(name,value) {
	if (!config.options[name])
		config.options[name] = value;
};

//
// Initializations
//

// require PasswordOptionPlugin 1.0.1 or better
bidix.checkPlugin("PasswordOptionPlugin", 1, 0, 1);

// styleSheet
setStylesheet('.txtUploadStoreUrl, .txtUploadBackupDir, .txtUploadDir {width: 22em;}',"uploadPluginStyles");

//optionsDesc
merge(config.optionsDesc,{
	txtUploadStoreUrl: "Url of the UploadService script (default: store.php)",
	txtUploadFilename: "Filename of the uploaded file (default: in index.html)",
	txtUploadDir: "Relative Directory where to store the file (default: . (downloadService directory))",
	txtUploadBackupDir: "Relative Directory where to backup the file. If empty no backup. (default: ''(empty))",
	txtUploadUserName: "Upload Username",
	pasUploadPassword: "Upload Password",
	chkUploadLog: "do Logging in UploadLog (default: true)",
	txtUploadLogMaxLine: "Maximum of lines in UploadLog (default: 10)"
});

// Options Initializations
bidix.initOption('txtUploadStoreUrl','');
bidix.initOption('txtUploadFilename','');
bidix.initOption('txtUploadDir','');
bidix.initOption('txtUploadBackupDir','');
bidix.initOption('txtUploadUserName','');
bidix.initOption('pasUploadPassword','');
bidix.initOption('chkUploadLog',true);
bidix.initOption('txtUploadLogMaxLine','10');


// Backstage
merge(config.tasks,{
	uploadOptions: {text: "upload", tooltip: "Change UploadOptions and Upload", content: '<<uploadOptions>>'}
});
config.backstageTasks.push("uploadOptions");


//}}}