Chapter 2:

numberOfBins = 100;
numberOfSteps = 100;
growthConstant = 3.6;
startValue = 0.5;
histogramR = zeros( 1, numberOfBins );
mapValuesR = zeros( 1, numberOfSteps );
mapValuesR( 1 ) = startValue;
histogramR( mapValuesR( 1 ) * numberOfBins + 1 ) = ...
      histogramR( mapValuesR( 1 ) * numberOfBins + 1 ) + 1;
for loop = 2 : numberOfSteps
	mapValuesR( loop ) = growthConstant * mapValuesR ...
      ( loop - 1 ) * ( 1 - mapValuesR( loop - 1 ) );
	histogramR( ceil( mapValuesR( loop ) * numberOfBins ) ) = ... 
      histogramR( ceil ( mapValuesR( loop ) * numberOfBins ) ) + 1;
end
figure( 1 )
plot( mapValuesR )
figure( 2 )
plot( histogramR ); 

Chapter 3:

#include <iostream.h>
#include "dislin.h"

const double KM = 1000;							
const double GRAVITATIONALCONSTANT = 6.67e-11;  
const double EARTHMASS = 5.97e24;
const double EARTHRADIUS = 6380 * KM;
const int MATSIZE = 20;			

double gravitationalField( double aX, double aY ) {
	double radius = sqrt( aX * aX + aY * aY )
	if ( radius <= EARTHRADIUS ) 				 
		return GRAVITATIONALCONSTANT * EARTHMASS * radius / 
			(  EARTHRADIUS * EARTHRADIUS * EARTHRADIUS ) ;
	else
		return GRAVITATIONALCONSTANT * EARTHMASS /
 			 ( aX * aX + aY * aY );
}

main ( ) {
	double position[matSize];		
	float field[matSize][matSize];
	float offset = matSize / 2 - 0.5; 
	for ( int loop = 0; loop < matSize; loop++ ) {
		position[loop] = 0.1 * earthRadius * (loop - offset);
	}
	float x, y;
	for ( int outerLoop = 0; outerLoop < matSize; outerLoop++ ) {
		x = position[outerLoop];
		for ( int innerLoop = 0; innerLoop < matSize; innerLoop++ ) {
			y = position[innerLoop];
			field[outerLoop][innerLoop] = gravitationalField( x, y );
		}
	}
	metafl( "XWIN" );
	disini( );		
	int iPlot = 2;		
	if ( iPlot == 1 )	
		qplsur( (float*) field, matSize, matSize );  
	else if (iPlot == 2)	
		qplclr( (float*) field, matSize, matSize );
	else {			
		int numberOfContours = 30;		
		qplcon( (float*) field, matSize, matSize, numberOfContours );
	}
}


Chapter 10

#include <strstream>
#include <iostream>
using namespace std;

class C {
      public:
	int iC;
	C( int aC = 0 ) : iC( aC ) { }
	string toString( ) {
		strstream s;
		s << "iC = " << iC <<ends;
		return s.str( );
	}
};

class D {
      public:
	C iC[2];
	int iD;
	D( int aD, C aC[ ] ) : iD( aD ) {
		iC[0] = aC[0];
		iC[1] = aC[1];
	}
	string toString( ) {
		strstream s;
		for ( int loop = 0; loop < 2; loop++ )
			s << iC[loop].toString( ) << " ";
		s << endl;
		s << "iD = " << iD;
		s << ends;
		return s.str( );
	}
};

main( ) {
	C C1[2] = {1, 2};
	D D1( 2, C1 );
	cout << D1.toString( );
}

Chapter 14

#include <sstream>

template <class C1, class C2 = int, int aN = 1> class Logical {
	C1 iC1;
	C2 iC2;
	int iN; 
	public:
		Logical ( C1 aC1, C2 aC2, int aN ) : iN( aN ) { iC1 = aC1, iC2 = aC2; }
		void and( ) { cout << ( iC1 && iC2 ) << " " << iN << endl; }
		string asString( ); 
};

template <class C1, class C2, int aN> string Logical<C1, C2, aN>::asString( ) {
	stringstream sout;
	sout << iC1 << ' ' << iC2 << endl;
	return sout.str( );
}

Chapter 15

#include <vector>
#include <iterator>

main ( ) {
	vector < int > aV( 3 );
	vector < int > :: iterator pItV;
	aV.push_back( 3 );	
	aV[1] = 1;
	aV.at( 2 ) =  2;		
	for ( pItV = aV.begin( ); pItV < aV.end( ); pItV++ )
		cout << *pItV << endl;
}


#include <set>						
#include <iterator>

main ( ) {
	double a[3] = { 1.3, 2.5, 0.3 };
	set < double > S1( a, a + 3 );
	set < double > :: iterator pItS;
	S1.insert( 0 );
	S1.insert( 1.3 );
	for ( pItS = S1.begin( ); pItS != S1.end( ); pItS++ )
		cout << *pItS << endl;
}


#include <set>
#include <iterator>

class C { 
	public:
	int value1, value2;
};

class myGreater {
	public:
	myGreater( ) {
	}
	bool operator( ) ( const C &C1, const C &C2 ) {
		return( C1.value1 < C2.value1 );
	}
};	

main ( ) {
	C C1 = { 1, 2 };
	C C2 = { 3, 4 };
	set <C, myGreater> s;
	set <C, myGreater> :: iterator pItV;
	s.insert( C1 );
	s.insert( C2 );
	for ( pItV = s.begin( ); pItV != s.end( ); pItV++ )
		cout << (*pItV).value1 << ' ';
}


#include <algorithm>
#include <iterator>
#include <vector>
#include <numeric> 	 					

int tripleFunction ( int aValue ) {return 3 * aValue;}

int tripleAccumulateFunction( int aPartialSum, int aValue ){
	return aPartialSum + 3 * aValue;
}

main( ) {
ostream_iterator <int> myOut( cout, " loop element \n" );
vector < int > myVector ( 20 );
fill( myVector.begin( ),  myVector.end( ), 5 ); 
replace( myVector.begin( ), myVector.begin( ) + 3, 5, 3 ); 
vector < int > tripleResult( 3 );
transform( myVector.begin( ) + 2, myVector.begin( ) + 4, tripleResult.begin( ),
 	tripleFunction ); 			 
sort( tripleResult.begin( ), tripleResult.end( ) );		
copy( tripleResult.begin( ), tripleResult.end( ), myOut );
cout << endl;
cout << accumulate( tripleResult.begin( ), tripleResult.end( ), 0 ); 
cout << endl;
cout << accumulate( tripleResult.begin( ), tripleResult.end( ), 0,    
      tripleAccumulateFunction );			
}

Chapter 16
 
<HTML> 
<TITLE> A Dislin Example </TITLE> 
<BODY> 
Text to appear before the applet 
<P><APPLET codebase="."code="dislin_1.class"width=900 height=700> 
Text to appear if an error prevents the applet from appearing 
</APPLET><P> 
Text to appear after the applet 
</BODY>  
</HTML>


import java.applet.Applet;
import java.awt.*;

public class myApplet extends Applet {
	private static Frame window = new Frame( "Drawing" );
	public void paint( Graphics g ) {
		g.drawString( "A Graph", 400, 400 );
		g.drawLine( 0, 10, 500, 500 );
		g.drawLine( 100, 100, 40, 600 );
		g.setColor( Color.cyan );
		g.fillOval( 300, 300, 80, 275 );
	}
}


Chapter 18

import java.awt.*;
import java.awt.event.*;

public class GraphExample extends Canvas implements MouseListener  {
    	public int iX = 0;
   	public int iY = 0;

    	private static Frame window = new Frame( "Drawing" );

    	public static void main ( String[ ] args ) {
		GraphExample g = new GraphExample( );
		g.setSize( 600, 800 );
		window.add( g );
		window.pack( );
		window.setVisible( true );
    	}

    	GraphExample( ) { addMouseListener( this );}

   	public void paint( Graphics g ) {
        		g.drawString( "A Graph", 400, 400 );
        		g.setColor( Color.cyan );
        		g.fillOval( 300, 300, 80, 275 );
        		g.drawLine( 100, 100, iX, iY );
    	}

    	public void mouseClicked( MouseEvent event ) {
		iX = event.getX( );
		iY = event.getY( );
		repaint( );
    	}
    	public void mouseExited( MouseEvent event ) {
		System.exit( 0 ); 
    	}

    	public void mouseEntered( MouseEvent event ) { }
    	public void mouseReleased( MouseEvent event ) { }
    	public void mousePressed( MouseEvent event ) { }

}


import java.awt.event.*;
import java.applet.Applet;
import java.awt.*;

class Beep implements ActionListener {
    	public void actionPerformed ( ActionEvent Event ) {
		Component C1 = (Component)Event.getSource( );
		C1.getToolkit( ).beep( );
		GraphExample2 g = new GraphExample2( );
    	}
}

public class GraphExample2 extends Applet  {
    	Point myPoint;
    	public void init( ) {
		Button MyButton = new Button( "Beep" );
		add( MyButton );
		MyButton.addActionListener( new Beep( ) );
    	}
    	public static Frame window = new Frame( "Drawing" );
    	public void paint( Graphics g ) {
		addMouseListener( new MouseAdapter( ) {    
	    		public void mouseClicked( MouseEvent event ) {
				myPoint = event.getPoint( );
				repaint( );
	    		}
	    		public void mouseExited( MouseEvent event ) {
				System.exit(0); 
	    		}
		} );
		g.drawString( "A Graph", 400, 400 );
		g.setColor( Color.cyan );
		g.fillOval( 300, 400, 80, 275 );
		g.drawLine( 100, 100, 300, 300 );
		g.drawLine( 100, 100, (int) myPoint.getX( ), (int) myPoint.getY( ) );
    	}
    	public void stop( ) {
		System.exit( 0 );
    	}
}   


class Thread1 extends Thread {
   	public void run( ) {
  		for (int loop = 0; loop < 10; loop++) {
         			System.out.println( loop );
		 	yield( );
		}
      }
}

class Thread2 extends Thread {
   	public void run( ) {
      		for (int loop = 0; loop < 10; loop++) {
         			System.out.println( loop );
		 	yield( );
		}
 	}
}

class TestThread {
   	public static void main( String args[ ] ) {
      		new Thread1( ).start( );
      		new Thread2( ).start( ); 
      }
}	


class MyStaticObject {
    	static Object O = new Object( );
}

class Thread1 extends Thread {
   	public  void run( ) {
   		synchronized ( MyStaticObject.O ) {
      			for (int loop = 0; loop < 10; loop++) {
         				System.out.println( loop );
		 		yield( );
			}
 		}
	}
}

class Thread2 extends Thread {
   	public  void run( ) {
   		synchronized ( MyStaticObject.O ) {
      			for (int loop = 0; loop < 10; loop++) {
         				System.out.println( loop );
		 		yield( );
			}
 		}
      }
}

class TestThread {
   	public static void main( String args[ ] ) {
      		new Thread1( ).start( );
      		new Thread2( ).start( ); 
      }
}


import java.io.*;
class MyClass implements Serializable {
   	private double[ ] myArray = { 0, 1, 2 };
   	public double getArrayElement( int aI ) { return 2 * myArray[aI]; } 
}
public class SerializableExample {
	public static void main( String args[ ] ) throws IOException {
   		MyClass MC1 = new MyClass( );
   		ObjectOutputStream OOS1 = new ObjectOutputStream( new 
			FileOutputStream( "myfile.ser" ) );
   		OOS1.writeObject( MC1 );
   		OOS1.close( ); 
      } 
}

import java.io.*; 
public class ReadSerializableExample {
      public static void main( String args[ ] ) throws IOException, 			
			ClassNotFoundException { 
		ObjectInputStream OI1 = new ObjectInputStream( new 		
			FileInputStream( "myfile.ser" ) ); 
   		MyClass MyObject = (MyClass) OI1.readObject( ); 
   		OI1.close( ); 
   		System.out.println( MyObject.getArrayElement( 0 ) ); 	
      } 
} 									


class MyObject<T> {
    	T iT;
    	void setT( T aT ) { iT = aT; }
}

class Generic {
      <C> void ClassType (C aC) { 
    		System.out.println( aC.getClass( ).getName( ) );
      }
      public static void main( String args[ ] ) {
    		MyObject<Double> MO1 = new MyObject<Double>( );
    		MO1.setT(20.0);
    		System.out.println( MO1.iT );		
		ClassType( new Double( 1.0 ) );		
      }
}

Chapter 19

double cube( double aD ) { 
      return pow( aD, 3 ); 
}

double linear( double aD ) { 
      return pow( aD, 1 ); 
} 

double derivOperator( double aF( double ), double aXValue, double aDeltaX ) {
return ( aF( aXValue + aDeltaX ) - aF( aXValue ) ) / aDeltaX; 
} 

main ( ) {
double deltaX = 1.e-1;
double xValue = 1.0;
int choice;
double (*myFunction) (double);
cout << "Choose a function  1 - cube, 2 - square  "<< endl;
cin >> choice;
switch ( choice ) {
      case 1: myFunction = cube;  break; 	 		
      case 2: myFunction = linear; break;  			
      default: cout << "Incorrect Input - program exiting"; 
}
      cout << derivOperator( myFunction, xValue , deltaX );  
}


main ( ) {
	double deltaX = 1.0e-1;
	double xValue = 1.0;
	float x[10], y[10]; 						
	for ( int loop = 0; loop < 10; loop++ ) {
		y[loop] = derivOperator( cube, xValue, deltaX );	
		x[loop] = deltaX;
		deltaX /= 2; 
	}
	metafl( "XWIN" );  	
	disini( );       					     
	name( "Step Length", "x" );				     	       
	name( "Error", "y" );  					     	       
	labels( "EXP","xy" );					     	       
	incmrk( 1 );						     	       
	setscl( x, 10, "X" );					     	       
	setscl( y, 10, "Y" );					     
	//scale( "LOG","XY" );					     
	float minX,maxX,minY,maxY,stepX,stepY;
	graf( minX, maxX,  minX, stepX, minY, maxY, minY, stepY ); 
	curve( x, y, 10 );            		
	disfin( );     				
} 


class DerivativeCalculator {
public: 
void setDx( double aDx ) { iDx = aDx; }
void setX( double aX ) { iX = aX; } 
double dx(  ) const { return iDx; } 
double calculateDerivative( ) const { return ( iF( iX + iDx ) - iF( iX )) / iDx; } 
DerivativeCalculator( double aX, double aDx, double aF( double ) ) : 
      iX( aX ), iDx( aDx ) { 
      iF = aF; 
      }
 	private:
double (*iF) ( double );    	 
double iDx;
      double iX; 
};

main ( ) {
double deltaX = 1.e-1;
double xValue = 1.0;
int choice;
double (*myFunction) ( double );
cout << "Choose a function  1 - cube, 2 - square  "<< endl;
cin >> choice;
switch ( choice ) {
      case 1: myFunction = cube;  break; 	 		
      case 2: myFunction = linear; break;  			
      default: cout << "Incorrect Input - program exiting"; 
}
      DerivativeCalculator DC1( 1.0,0.1, myFunction );
      cout << DC1.calculateDerivative( );
}

function result = myMidpoint( aFunction, leftEndPoint, rightEndPoint, numberOfIntervals )
deltaX = ( rightEndPoint - leftEndPoint ) / numberOfIntervals;
result = deltaX * sum( aFunction( leftEndPoint + 0.5 * deltaX + ...
       		deltaX * ( 0 : numberOfIntervals - 1 ) ) )

 
function result = mySimpson( aFunction, aLeftEndPoint, aRightEndPoint, ... 
	aNumberOfIntervals )
deltaX = ( aRightEndPoint - aLeftEndPoint ) / aNumberOfIntervals;
xLeft = leftEndPoint
deltaX2 = deltaX / 2;
mySum = 0;
for loop = 0 : aNumberOfIntervals - 1
    	mySum = mySum + 2 * aFunction( xLeft  ) + ...
      4 * aFunction( xLeft + deltaX2 );
    	xLeft = xLeft + deltaX;
end
mySum = mySum + aFunction( aRightEndPoint ) - aFunction( aLeftEndPoint );
result = deltaX * mySum / 6; 


function result = myBisect( aFunction, aLeftLimit, aRightLimit, aError )
while ( abs( aLeftLimit - aRightLimit ) > aError ) 
    	midpoint = ( aLeftLimit + aRightLimit ) / 2;
    	if aFunction ( midpoint ) * aFunction( aRightLimit ) <= 0
        		aLeftLimit = midpoint;
    	else
        		aRightLimit = midpoint;
    	end
end
result = ( aLeftLimit + aRightLimit ) / 2; 


interface MyFunction {
    	double  myFunction( double aD );
}

public class Newton {
    	static int count;
    	static double derivative( MyFunction aMF, double aX ) {
		double deltaX = 1.e-3;
		return ( aMF.myFunction( aX + deltaX ) 
			- aMF.myFunction( aX - deltaX ) ) / ( 2. * deltaX ); 
    	}
    	static double newton( MyFunction aMF1, double aEstimate, double aError ){
		double deltaX = - aMF1.myFunction( aEstimate ) 
			/ derivative( aMF1, aEstimate );
		aEstimate += deltaX;
		if ( Math.abs( deltaX ) < aError || count++ == 40 ) return aEstimate; 
		else return newton( aMF1, aEstimate, aError ); 
    	}
    	public static void main  ( String [ ] args ) throws Exception {
		System.out.println( "Insert 1 for square, 2 for cube" );
		char c = (char) System.in.read( );
		MyFunction MF1;
		switch ( c ) {
	    		case '1': MF1 = new MyFunction( ) { 
				public double myFunction( double aD ) 
				{ return aD * aD - 4; } } ; 
			break;
	    		case '2': MF1 = new MyFunction( ) 
				{ public double myFunction( double aD ) 
				{ return aD * aD * aD - 27; } } ; 
			break;
	    		default: System.out.println( "Error in input" ); return;
		}
		double estimate = 2.5;
		double error = 1.e-8;
		System.out.println( newton( MF1, estimate, error ) );
    	}
}


function result = myMinimum( aFunction, aLeftEndPoint, aRightEndPoint, ...
      aError );
leftValue = aFunction( aLeftEndPoint );
rightValue = aFunction( aRightEndPoint );
for loop = 1 : 500
	middlePoint = ( aLeftEndPoint + aRightEndPoint ) / 2;
	middleValue = aFunction( middlePoint );
	aNewPoint = ( aRightEndPoint + middlePoint ) / 2;
	newValue = aFunction( aNewPoint );
	if ( newValue < middleValue )
		aLeftEndPoint = middlePoint;
	else
		aRightEndPoint = aNewPoint;
	end
	aNewPoint = ( aLeftEndPoint + middlePoint ) / 2;
	newValue = aFunction( aNewPoint );
	if ( newValue < middleValue )
		aRightEndPoint = middlePoint;
	else
		aLeftEndPoint = aNewPoint;
	end
	if aRightEndPoint - aLeftEndPoint < aError
		break;
	end
end
result = middlePoint;

 
#include <stdio.h>
const int n = 2;
void gauss( double aA[ ][n], double aC[ ], double aB[  ] ) {
	for ( int i = 0; i < n; i ++ ) {
		if ( !aA[i][i] ) exit( 0 );
		for ( int j = i + 1; j < n; j++ ) {
			double d = aA[j][i] / aA[i][i];
			for ( int k = i + 1; k < n; k++ ) aA[j][k] -= d*aA[i][k];
			aB[j] -= d*aB[i];
		}
	}
	if ( !aA[n-1][n-1] ) exit( 0 );
	for ( int i = n - 1; i >= 0; i-- ) {
		aC[i] = aB[i];
		for ( int j = i + 1; j < n; j++ )
			aC[i] -= aA[i][j] * aC[j];
		aC[i] /= aA[i][i];
	}
}

main( ) {
	double a[n][n] = { { 1, 2 }, { 3, 8 } };
	double b[n] = { 2, 5 };
	double c[n];
	gauss( a, c, b );
	cout << c[0] << '\t' << c[1] << endl;
}


void tridiagonalSolver( double aLowerCodiagonal[ ], double aDiagonal[ ], 
                       double aUpperCodiagonal[ ], double aInputVector[ ],
                       double aOutputVector[ ], double aScratch[ ], 
                       int aNumberOfPoints ) {
    	double bsave = aDiagonal[0];
    	if ( !bsave ) exit( 0 );
    	aOutputVector[0] = aInputVector[0] / bsave;
    	for ( int loop = 1; loop < aNumberOfPoints; loop++ ) {
        		aScratch[loop] = aUpperCodiagonal[loop - 1] / bsave;
        		bsave = aDiagonal[loop] - aScratch[loop] * aLowerCodiagonal[loop - 1];
        	if ( !bsave ) exit( 0 );
        	aOutputVector[loop] = ( aInputVector[loop] - 
                              	aLowerCodiagonal[loop - 1] * aOutputVector[loop - 1] ) 
                              	/ bsave;
    	}
    	for ( int loop = aNumberOfPoints - 2; loop > -1; loop-- )
        		aOutputVector[loop] - = aScratch[loop + 1] 
      * aOutputVector[loop + 1];
}

main( ) {
    double diagonal[2] = { 1, 8 };
    double upperCodiagonal[1] = { 2 };
    double lowerCodiagonal[1] = { 3 };
    double inputVector[2] = { 2, 5 };
    double outputVector[2];
    double scratch[2];
    tridiagonalSolver( lowerCodiagonal, diagonal, upperCodiagonal, 
                      inputVector, outputVector, scratch, 2 );
    cout << outputVector[0] << '\t' << outputVector[1] << endl;
}


function outputVector = myTridigonal( aLowerCodiagonal, ...
      aDiagonal, aUpperCodiagonal, aInputVector )

numberOfEquations = length( aDiagonal );
outputVector = zeros(1, numberOfEquations);
for loop = 2 : numberOfEquations 
        temporary = aLowerCodiagonal(loop - 1) / aDiagonal(loop - 1);
        aDiagonal(loop) = aDiagonal(loop) - temporary * aUpperCodiagonal(loop - 1);
        aInput(loop) = aInputVector(loop) - temporary * aInput(loop - 1);
end
outputVector(numberOfEquations) = aInput(numberOfEquations) ...
       / aDiagonal(numberOfEquations); 
for loop = numberOfEquations - 1 : -1 : 1
        outputVector(loop) = ( aInput(loop) - aUpperCodiagonal(loop) ...
      * outputVector(loop + 1) ) / aDiagonal(loop); 
end


const NUMBEROFTIMESTEPS = 100;
main( ) {
      double x[NUMBEROFTIMESTEPS], v[NUMBEROFTIMESTEPS],
 		 k=1, m=1, dt=0.06;
      x[0] = 0;
      v[0] = 1; 			
      for ( int loop = 1; loop < NUMBEROFTIMESTEPS; loop++ ) {
	x[loop] = x[loop - 1] + dt * v[loop - 1]; 
	v[loop] = v[loop - 1] - k * dt * x[loop - 1] / m; 
      } 
      qplot( x, v, NUMBEROFTIMESTEPS )
}

 
numberOfTimeSteps = 100;
deltaTime = 0.2;
gravitationalConstant = -9.8;
dragConstant = 6.0E-2;

Ball.positionRC = zeros( 2, numberOfTimeSteps );
Ball.velocityRC = zeros( 2, numberOfTimeSteps ); 
Ball.velocityRC(:, 1) = [ 10; 100 ];		

for loop = 2 : numberOfTimeSteps
	force = [ 0 ; gravitationalConstant ] - + ...
		dragConstant * Ball.velocityRC(:, loop - 1); 
	Ball.positionRC(:, loop) = Ball.positionRC(:, loop - 1) + ...
		deltaTime * Ball.velocityRC(:, loop - 1);
	Ball.velocityRC(:, loop) = Ball.velocityRC(:, loop - 1) + ...
		deltaTime * force;
end

plot( Ball.positionRC(1, :), Ball.positionRC(2, :));


f = @cos;
numberOfPoints = 10000;
lowerLimit = 0;                
upperLimit = pi / 2;           
yLower = -2;          
yUpper = +2;          
regionArea = abs( yUpper - yLower ) * ( upperLimit - lowerLimit );
xValuesR = rand( 1, numberOfPoints ) * ( upperLimit - lowerLimit ) + lowerLimit;
yValuesR = rand( 1, numberOfPoints ) * ( yUpper - yLower ) + yLower;
numberOfPointsBelowCurve = sum( yValuesR < f( xValuesR ));
integral = numberOfPointsBelowCurve * regionArea / numberOfPoints + ...
     yLower * ( upperLimit - lowerLimit )


clear all
numberOfSteps = 40;
numberOfRealizations = 50000;
histogramR = zeros( 1, numberOfSteps + 1 );
for loop = 1 : numberOfRealizations
    	histogramIndex = sum( round( rand( 1, numberOfSteps ) ) ) + 1;
    	histogramR( histogramIndex ) = histogramR( histogramIndex ) + 1;
end
histogramR = histogramR / sum( histogramR );
xScaleR = 2 * ( 0 : numberOfSteps ) - numberOfSteps;
semilogy( xScaleR, histogramR, "o", 'markersize', 3 );

clear all
numberOfSteps = 40;
numberOfRealizations = 50000;
histogramR = zeros( 1, numberOfSteps + 1 );
for loop = 1 : numberOfRealizations
    	histogramIndex = sum( round( rand( 1, numberOfSteps ) ) ) + 1;
    	histogramR( histogramIndex ) = histogramR( histogramIndex ) + 1;
end
for loop = 1 : numberOfSteps;
    	histogramR( loop ) = histogramR( loop ) * ...
		0.4^loop * 0.6^( numberOfSteps - loop );
end
histogramR = histogramR / sum( histogramR );
xScaleR = 2 * ( 0 : numberOfSteps ) - numberOfSteps;
semilogy( xScaleR, histogramR, 'x', 'markersize', 3 );

clear all
numberOfRealizations = 100000;
numberOfSteps = 40;				
myBeta = 1.0;
histogramR = zeros( 1, numberOfSteps + 1 );
stepSequenceR = round( rand( 1, numberOfSteps ) );
histogramIndex = sum( stepSequenceR ) + 1;
centralPoint = numberOfSteps / 2;
for loop = 1 : numberOfRealizations;
    	flipPosition = fix( rand * numberOfSteps ) + 1;
    	stepSequenceR(flipPosition) = 1 - stepSequenceR(flipPosition);
    	histogramIndexNew = sum( stepSequenceR ) + 1;
    	if rand < exp( myBeta * ( abs( histogramIndexNew - centralPoint ) - ...
            	abs( histogramIndex - centralPoint ) ) )
       		histogramIndex = histogramIndexNew;
    	else
       		stepSequenceR(flipPosition) = 1 - stepSequenceR(flipPosition);
    	end  
    	histogramR(histogramIndex) = histogramR(histogramIndex) + 1;
end
histogramR = histogramR .* exp( - myBeta * ... 
      abs( [1 : numberOfSteps + 1] - centralPoint ) );
xScaleR = 2 * ( 0 : numberOfSteps ) - numberOfSteps;
semilogy( xScaleR, histogramR / sum( histogramR ) );


clear all
numberOfRealizations = 100000;
numberOfSteps = 40;
numberOfOuterLoops = 4;
histogramR = ones( 1, numberOfSteps + 1 );
stepSequenceOldR = round ( rand( 1, numberOfSteps ) );
stepSequenceNewR = stepSequenceOldR;
histOld = 1;
histogramIndexOld = sum( stepSequenceOldR ) + 1;
 
for loopOuter = 1 : numberOfOuterLoops; 
    	histogramNewR = ones( 1, numberOfSteps + 1 );
    	for loop = 1 : numberOfRealizations;       
        	stepSequenceNewR = stepSequenceOldR;
        	flipPosition = ceil( numberOfSteps * rand );
        	stepSequenceNewR(flipPosition) = ...
			1 - stepSequenceOldR( flipPosition );
        	histogramIndexNew = sum( stepSequenceNewR ) + 1;
        	histNew = histogramR( histogramIndexNew );
        	if ( rand < histOld / histNew )	
			stepSequenceOldR = stepSequenceNewR;
            		histOld = histNew;
			histogramIndexOld = histogramIndexNew;
        	end      
        	histogramNewR( histogramIndexOld ) = ...
			histogramNewR( histogramIndexOld ) + 1;        
    	end
    	histogramR = histogramR .* histogramNewR;
end
 
xScaleR = 2 * ( 0 : numberOfSteps ) - numberOfSteps;
semilogy( xScaleR, histogramR / sum( histogramR ) );

clear all
numberOfParticles = 20;
numberOfCollisions = 200;
windowSize = 500;
collisionTime = 0;
Particle.positionR = rand( 1, numberOfParticles ) * windowSize;
Particle.positionR( : ) = sort( Particle.positionR( : ) );
Particle.velocityR = ( rand( 1, numberOfParticles ) - 0.5 ) * 2;
for loop = 1 : numberOfCollisions
    minimumTime = 1.e20;
    for loopParticle = 2 : numberOfParticles
        if Particle.velocityR(loopParticle - 1) - Particle.velocityR(loopParticle) > 0
            deltaTime = ( Particle.positionR(loopParticle) - 
		Particle.positionR(loopParticle - 1) ) / ...
                ( Particle.velocityR(loopParticle - 1) - Particle.velocityR(loopParticle) );
            if  deltaTime < minimumTime
                minimumTime = deltaTime; 
                particleIndex = loopParticle;
                isWall = 0;
            end;
        end
    end
    if Particle.velocityR(numberOfParticles) > 0
        deltaTime = ( windowSize - Particle.positionR(numberOfParticles) ) / ...
            Particle.velocityR(numberOfParticles);
        if  deltaTime < minimumTime; 
            minimumTime = deltaTime; 
            particleIndex = numberOfParticles;
            isWall = 1;
        end;
    end
    if Particle.velocityR(1) < 0
        deltaTime = - Particle.positionR(1) / Particle.velocityR(1);
        if  deltaTime < minimumTime; 
            minimumTime = deltaTime; 
            particleIndex = 1;
            isWall = 1;
        end;
    end;
    for loopParticle = 1 : numberOfParticles
        Particle.positionR(loopParticle) = Particle.positionR(loopParticle) + ...
            Particle.velocityR(loopParticle) * minimumTime;
        if loopParticle == particleIndex
            if isWall 
                Particle.velocityR(loopParticle) = -Particle.velocityR(loopParticle);
            else
                saveVelocity = Particle.velocityR(loopParticle);
                Particle.velocityR(loopParticle) = Particle.velocityR(loopParticle - 1);
                Particle.velocityR(loopParticle - 1) = saveVelocity;
            end
        end
    end
    collisionTime = collisionTime + minimumTime;
    saveTimeR(loop) = collisionTime;
    savePositionR(loop) = Particle.positionR(numberOfParticles / 2);
end
plot( saveTimeR, savePositionR );


clear all
Ising.numberOfEnergies = 1;
Ising.energy(1) = 0;
Ising.numberOfStates(1) = 0;
numberOfSpins = 8;
numberOfRealizations = 2^numberOfSpins;
for loop = 1 : numberOfRealizations
	spins = dec2binase( loop - 1, 2, numberOfSpins );
	energy = 0;
	for innerLoop = 1 : numberOfSpins
            energy = energy - str2num( spins(innerLoop) ) * ...
                 str2num( spins( mod( innerLoop, numberOfSpins ) + 1 ) );
	end
	energyFound = 0;
	for innerLoop = 1 : Ising.numberOfEnergies
		if Ising.energy(innerLoop) == energy;
			Ising.numberOfStates(innerLoop) = ... 
	 			Ising.numberOfStates(innerLoop) + 1;
			energyFound = 1;
		end
		if energyFound == 1; break; end;
		energyFound
	end
	if energyFound == 0
		Ising.numberOfEnergies = Ising.numberOfEnergies + 1;
		Ising.energy(Ising.numberOfEnergies) = energy;
		Ising.numberOfStates(Ising.numberOfEnergies) = 1;
	end
end
myBeta = 0.05;
partitionFunction = 0;
for loop = 1 : Ising.numberOfEnergies
	partitionFunction = partitionFunction + Ising.numberOfStates(loop) * ...
		exp( -myBeta * Ising.energy(loop) );
end
averageEnergy = 0;
averageSquaredEnergy = 0;
for loop = 1 : Ising.numberOfEnergies
	averageEnergy = averageEnergy + Ising.energy(loop) * ...
		Ising.numberOfStates(loop) * ...
		exp( -myBeta * Ising.energy(loop) );
	averageSquaredEnergy = averageSquaredEnergy + ...
		Ising.energy(loop)^2 * Ising.numberOfStates(loop) * ...
		exp( -myBeta * Ising.energy(loop) );
end
averageEnergy = averageEnergy / partitionFunction;
averageSquaredEnergy = averageSquaredEnergy / partitionFunction;
specificHeatPerSpin = myBeta^2 * ( averageSquaredEnergy - averageEnergy^2 ) / numberOfSpins; 

clear all
numberOfPoints = 6;			
halfWidth = 5;				
xPositionR = linspace( -halfWidth, halfWidth, numberOfPoints );
yPositionR = xPositionR;
field.xValueRC = zeros( numberOfPoints, numberOfPoints );
field.yValueRC = zeros( numberOfPoints, numberOfPoints );

for outerLoop = 1 : numberOfPoints
	for innerLoop = 1 : numberOfPoints
		radius = sqrt( xPositionR(outerLoop)^2 + yPositionR(innerLoop)^2 );
		unitX = xPositionR(outerLoop) / radius;
		unitY = yPositionR(innerLoop) / radius;
		field.xValueRC(innerLoop, outerLoop) = unitX / radius^2;
		field.yValueRC(innerLoop, outerLoop) = unitY / radius^2;
	end
end

quiver( xPositionR, yPositionR, field.xValueRC, field.yValueRC ); 
 

const int NUMBEROFPOINTS = 100;

main( ) {
	double deltaTime = 0.01;
	int numberOfTimeSteps = 5000;
	double deltaX = 0.1;
	double diffusionConstant = 0.5;
	double coefficient  = diffusionConstant * deltaTime / ( deltaX * deltaX );

	float position[NUMBEROFPOINTS];
	float tempLast[NUMBEROFPOINTS], temp[NUMBEROFPOINTS ] = { 0 };
	for ( int loop = 0; loop < NUMBEROFPOINTS; loop++ ) {
		position[ loop ] = deltaX * loop;
		temp[loop] = sin( loop * M_PI / ( NUMBEROFPOINTS - 1 ) );
	}

	for ( int outerLoop = 0; outerLoop < numberOfTimeSteps; outerLoop++ ) {
		if ( outerLoop % 500 == 0 ) qplot( position, temp, NUMBEROFPOINTS );
		for ( int loop = 1; loop < NUMBEROFPOINTS - 1; loop++ )
			tempLast[loop] = temp[loop];
		for ( int loop = 1; loop < NUMBEROFPOINTS - 1; loop++ ) 
			temp[loop] = coefficient * ( tempLast[loop - 1] - 2 * 
			tempLast[loop] + tempLast[loop + 1] ) + temp[loop];
	}
} 


import de.dislin.*;

public class Hopscotch {

    static public void main( String[ ] args ) {
	int numberOfPoints = 200;
	int numberOfSteps = 500;
	float deltaX = 1.0F;
	float velocity = 1.0F;
	float deltaT = 1.0F;
	float field1[ ] = new float[numberOfPoints];
	float field2[ ] = new float[numberOfPoints];
	float xValues[ ] = new float[numberOfPoints];
	for ( int loop = 0; loop < numberOfPoints; loop++ ) {
	    xValues[loop] = ( (float) loop ) * deltaX;
	}
	for ( int outerLoop = 0; outerLoop < numberOfSteps - 1; outerLoop++ ) { 
	    for ( int loop = 0; loop < numberOfPoints - 1; loop++ ) 
		field2[loop] += velocity * deltaT * ( field1[loop + 1] - field1[loop] ) / deltaX ;
	    for ( int loop = 1; loop < numberOfPoints; loop++ ) 
		field1[loop] += velocity * deltaT * ( field2[loop] - field2[loop - 1] ) / deltaX ;
	    field1[0] = (float) Math.exp( - ( outerLoop - 50 ) * ( outerLoop - 50 ) / 80. );
	    if( outerLoop % 10 == 0 ) Dislin.qplot( xValues, field1, numberOfPoints );
	}
      }
}


vacuumWaveLength = 1.0;
gridLength = 200;
numberOfGridPoints = 1024;		
numberOfSpatialSteps = 20;		
propagationDistance = 1000;		
referenceIndex = 1.0;
propagationStepLength = propagationDistance / numberOfSpatialSteps;
pointLocationsPlus1R = linspace( -gridLength / 2, gridLength / 2, ... 
      numberOfGridPoints + 1);
pointLocationsR = pointLocationsPlus1R( 1 : numberOfGridPoints );
k0 = 2 * pi / vacuumWaveLength;		
fieldR = exp( -pointLocationsR.^2 / ( 2. * 5.0^2 ) );
squaredFourierWavevectorComponentsR =  - ( 2 * pi / gridLength )^2 * ...
      [ 0 : numberOfGridPoints / 2, -numberOfGridPoints / 2 + 1 : -1 ].^2;
propagationOperatorR = exp( - i * propagationStepLength / ...
      ( 2 * k0 * referenceIndex ) * squaredFourierWavevectorComponentsR );
lensOperatorR = exp( i * k0 * pointLocationsR.^2 / 500 );
for loop = 1 : numberOfSpatialSteps
      fieldR = ifft( propagationOperatorR .* fft( fieldR ) );
      if loop == floor( numberOfSpatialSteps / 2 )
      fieldR = lensOperatorR .* fieldR;
      end
      plot( pointLocationsR, abs( fieldR ) )
      drawnow
end


gridLength = 50;
numberOfGridPoints = 512;			
numberOfTimeSteps = 500;
propagationTime = 8;
propagationStepTime = propagationTime / numberOfTimeSteps;
pointLocationsPlus1R = linspace( -gridLength / 2, gridLength / 2, ... 
      numberOfGridPoints + 1);
pointLocationsR = pointLocationsPlus1R( 1 : numberOfGridPoints );
fieldR = 1.0 ./ cosh( pointLocationsR - 8. ) .* ...
	exp( -i * 2 * pi * sin( 50 * pi / 180. ) * pointLocationsR );
propagationOperatorR = exp( -i * propagationStepTime * ( 2 * pi / gridLength )^2 / 4 * ...
		[ 0 : numberOfGridPoints / 2 , -numberOfGridPoints / 2 + 1 : -1 ] .^2 );
for loop = 1 : numberOfTimeSteps
	fieldR = ifft( propagationOperatorR .* ...
		fft( exp( i * propagationStepTime * conj( fieldR ) .* fieldR ) ...
		.* ifft( propagationOperatorR .* fft( fieldR ) ) ) );
	if rem( loop - 1, 50 ) == 0 
   		plot( pointLocationsR, abs( fieldR ) )
   		drawnow
	end
end


function y = potential ( x )
wellWidth = 5;
wellDepth = 5;
y = zeros( length( x ) , 1 );
for loop = 1 : length( x ) 
	if abs( x(loop) ) < wellWidth;
		y(loop) = -wellDepth;
	end
end

hold on;
stepLength = 0.005;
numberOfTimeSteps = 500;
fieldWidth = 2;
computationalWindowWidth = 20;
numberOfPoints = 300;			
gridPointsPlus1R = linspace( -computationalWindowWidth / 2, computationalWindowWidth / 2, numberOfPoints + 1);
gridPointsC = gridPointsPlus1R( 1 : numberOfPoints ).';
wavefunctionC = exp( -gridPointsC.^2 ./ ( 2 * fieldWidth^2 ) );
propagationOperatorC = exp( -i * stepLength / 2 * ...
	( 2 * pi / computationalWindowWidth )^2 * ...
	([ 0 : numberOfPoints / 2 , -numberOfPoints / 2 + 1 : -1 ]').^2 );
	potentialOperatorC = exp( -i * stepLength * potential( gridPointsC ) );
for loop = 1 : numberOfTimeSteps
  	wavefunctionC = ifft( propagationOperatorC .* fft( potentialOperatorC .* wavefunctionC ) );
    	if ( rem(loop,50) == 0 )
       		plot( gridPointsC, abs( wavefunctionC ), 'r' );
                        coefficient = ( 1 + ( stepLength * loop )^2 / fieldWidth^4 );
        		plot( gridPointsC, max( abs( wavefunctionC ) )coefficient^-0.25 * ...
				exp( -gridPointsC.^2 ./ ( 2 * ( fieldWidth^2 * coefficient ) ), ‘k’ ); 
    	end
    	drawnow
end

hold on;
stepLength = .005;
numberOfTimeSteps = 500;
fieldWidth = 2;
computationalWindowWidth = 20;
numberOfPoints = 300;
dx = computationalWindowWidth / ( numberOfPoints - 1 );
gridPointsC = linspace( -computationalWindowWidth / 2, ...
    	computationalWindowWidth / 2, numberOfPoints ).';
wavefunctionC = exp( -gridPointsC.^2 ./ ( 2 * fieldWidth^2 ) );
enableMyTridiagonal = 1;
if enableMyTridiagonal
    	aC = stepLength/( 4 * 1i ) * ( ones( numberOfPoints - 1, 1 ) ) / ( 2 * dx * dx );
    	bC = 0.5 * ones( numberOfPoints, 1 ) - stepLength/( 4 * 1i ) * ... 
		( potential( gridPointsC ) + ones( numberOfPoints, 1 ) / ( dx * dx ) );
else
    	onesVectorC = ones( numberOfPoints, 1 );
    	leftMatrixRCs = speye( numberOfPoints ) + i * stepLength / 2 * ...
        		( spdiags( potentialC, 0, numberOfPoints, numberOfPoints ) - ...
        		spdiags( [ onesVectorC, -2 * onesVectorC, onesVectorC ] / ...
        		( 2 * dx * dx ), -1 : 1, numberOfPoints, numberOfPoints ) );
end
for loop = 1 : numberOfTimeSteps
    	if enableMyTridiagonal
        		wavefunctionC = myTridiagonal( aC, bC, aC, wavefunctionC ) - wavefunctionC;
    	else
        		wavefunctionC = 2 * ( leftMatrixRCs \ wavefunctionC ) - wavefunctionC;
    	end
    	if ( rem(loop,50) == 0 )
       		plot( gridPointsC, abs( wavefunctionC ), 'g' );
    	end
    	drawnow
end


function outputVector = myTridiagonal( aLowerCodiagonal, aDiagonal, aUpperCodiagonal, aInputVector )
numberOfEquations = length( aDiagonal );
rowSize = size( aLowerCodiagonal );
if rowSize == 1
      outputVector = zeros(1, numberOfEquations)
else
      outputVector = zeros( numberOfEquations, 1 ); 
end
for loop = 2 : numberOfEquations
    	temporary = aLowerCodiagonal( loop - 1 ) / aDiagonal( loop - 1 );
    	aDiagonal( loop ) = aDiagonal( loop ) - temporary * ...
        	aUpperCodiagonal( loop - 1 );
    	aInputVector( loop ) = aInputVector( loop ) - temporary * ...
        	aInputVector( loop - 1 );
end
outputVector( numberOfEquations ) = aInputVector( numberOfEquations ) / ...
    	aDiagonal( numberOfEquations );
for loop = numberOfEquations - 1 : -1 : 1
    	outputVector( loop ) = ( aInputVector( loop ) - aUpperCodiagonal( loop ) * ...
        	outputVector( loop + 1 ) ) / aDiagonal( loop );
end

numberOfPoints = 100;
computationalWindowWidth = 20;
deltaX = computationalWindowWidth / ( numberOfPoints - 1 );
leftEndPoint = -computationalWindowWidth / 2;
MRC = zeros( numberOfPoints, numberOfPoints );
xR = linspace( leftEndPoint, -leftEndPoint, numberOfPoints );
for ( loop = 1 : numberOfPoints )
	MRC(loop, loop) = 1 / deltaX^2 + potential( xR(loop) );
	if ( loop != numberOfPoints ) 
		MRC(loop, loop + 1) = -1 / ( 2 * deltaX^2 ); 
	end;
	if ( loop != 1 ) 
		MRC(loop, loop - 1) = -1 / ( 2 * deltaX^2 ); 
	end;
end
[ eigenVectorsC, eigenValues ] = eigs( MRC, 1, 'sa' );
plot( xR, eigenVectorsC * sign( eigenVectorsC(numberOfPoints / 2) ) );

numberOfPoints = 100;
computationalWindowWidth = 20;
deltaX = computationalWindowWidth / ( numberOfPoints - 1 );
leftEndPoint = -computationalWindowWidth / 2;
SRC = zeros( numberOfPoints, numberOfPoints );
SBarRC = zeros( numberOfPoints, numberOfPoints );
xR = linspace( leftEndPoint, -leftEndPoint, numberOfPoints );
for ( loop = 1 : numberOfPoints )
	SRC(loop, loop) = 1 / ( 2 * deltaX ) + potential( xR(loop) ) * deltaX / 3;
	SBarRC(loop, loop) = deltaX / 3;
	if ( loop ~= numberOfPoints)
		SRC(loop, loop + 1) =  -1 / ( 2 * deltaX ) + potential( xR(loop) ) * deltaX / 6;
		SBarRC(loop, loop + 1) = deltaX / 6;
	end	
end
SRC = SRC + SRC';
SBarRC = SBarRC + SBarRC';
[ eigenVectorsC, eigenValues ] = eigs( SRC, SBarRC, 1, 'sa' );
plot( xR, eigenVectorsC * sign( eigenVectorsC(numberOfPoints / 2) ) );