added another example

doc/src/How2ReadData/How2ReadData.do.txt

@@ -1041,3 +1041,51 @@ p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}\left(c_{ij}+1\right)^{\gamma},
 !et
 where $c_{ij}$ represents the number of  previous interactions that have taken place between $i$ and $j$. The factor $1$ is added in order to ensure that if they have not interacted earlier they can still interact.  Perform similar studies as above with $N=1000$, $\alpha=1.0$ and $\alpha=2.0$ using $\gamma = 0.0, 1.0, 2.0, 3.0$ and $4.0$. Plot the wealth distributions for these cases and try to extract eventual power law tails with and without a saving $\lambda$ in each transaction. Comment your results and compare them with figures 5 and 6 of "Goswami and Sen":"http://www.sciencedirect.com/science/article/pii/S0378437114006967".  
 
+!split
+===== Particle in one dimension an velocity distribution =====
+!bc pycod 
+# Program to test the Metropolis algorithm with one particle at given temp in one dimension
+import numpy as np
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+import random
+from math import sqrt, exp, log
+# initialize the rng with a seed
+random.seed()
+# Hard coding of input parameters
+MCcycles = 100000
+Temperature = 2.0
+beta = 1./Temperature
+InitialVelocity = -2.0
+CurrentVelocity = InitialVelocity
+Energy = 0.5*InitialVelocity*InitialVelocity
+VelocityRange = 10*sqrt(Temperature)
+VelocityStep = 2*VelocityRange/10.
+AverageEnergy = Energy
+AverageEnergy2 = Energy*Energy
+VelocityValues = np.zeros(MCcycles)
+# The Monte Carlo sampling with Metropolis starts here
+for i in range (1, MCcycles, 1):
+    TrialVelocity = CurrentVelocity + (2.0*random.random() - 1.0)*VelocityStep
+    EnergyChange = 0.5*(TrialVelocity*TrialVelocity -CurrentVelocity*CurrentVelocity);
+    if random.random() <= exp(-beta*EnergyChange):
+        CurrentVelocity = TrialVelocity
+        Energy += EnergyChange
+        VelocityValues[i] = CurrentVelocity
+    AverageEnergy += Energy
+    AverageEnergy2 += Energy*Energy
+#Final averages
+AverageEnergy = AverageEnergy/MCcycles
+AverageEnergy2 = AverageEnergy2/MCcycles
+Variance = AverageEnergy2 - AverageEnergy*AverageEnergy
+print(AverageEnergy, Variance)
+n, bins, patches = plt.hist(VelocityValues, 400, facecolor='green')
+
+plt.xlabel('$v$')
+plt.ylabel('Velocity distribution P(v)')
+plt.title(r'Velocity histogram at $k_BT=2$')
+plt.axis([-5, 5, 0, 600])
+plt.grid(True)
+plt.show()
+
+!ec