AppendixB.tex

\chapter{Including Your Code} \label{AppendixCode}
\label{app:Talk}

You might want to include some code in your thesis. This is done with a lstlisting environment. Using lstlisting provides syntax highlighting. 

\begin{lstlisting}[language=Python] 
x=10
y=x**2
print('Helo world', y)
\end{lstlisting}

A whole code might look like this

\small
\begin{lstlisting}[language=Python]
#############################################################
# Program to do a curve fit in python with a user defined 
#  fit equation and %  with data supplied by the user.
# The user types the data into numpy arrays by hand,
#  and supplies a function with the fit equation.
#  In this example, the function is called RC_Charging.
#
# Author:  Todd Lines
# Date:    2023-02-20
#############################################################
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Define the Gaussian function
print("  ")
print("Find a curve fit to a user defined function")

#Here is where you define the function to use in the fit
def RC_Charging(x, A, B):
x=np.longdouble(x)
y = A*(1-np.exp(-B*x))
return y

#here is where you give the data to fit. Put it into numpyu arrays
xdata = np.array([10,    20,   30,   40,  50,   60,   70,   80,   
90,  100, 110, 120,  130,  140,  150,  160,  
170,  180])
ydata = np.array([0.5, 0.76, 0.97, 1.15, 1.3, 1.42, 1.54, 1.64, 
1.73, 1.82, 1.9, 1.97,2.05, 2.12, 2.18, 2.24, 
2.28, 2.33])

#Now plot the data so we can see the data points
plt.plot(xdata, ydata, 'o')


#Now perform the curve fit. We can't just use a linear fit. 
#  The data is very much #  not linear. So we will use a more 
#  robust curve fit rotine from scipy. The sciepy optimize 
#  curve_fit() routine needs the equation to use for the
#  fit as a function (here RC_Charting()) and it needs the 
#  fit parameters and for the error on 
#  the fit parameters we need the covariance matrix to be output
parameters, covariance = curve_fit(RC_Charging, xdata, ydata)

#Pull out the fit prameters
fit_A = parameters[0]
fit_B = parameters[1]


#Pull out the uncertanty from the diagonal elements of the 
#  covariance matrix. Remember that the diagonal elements 
#  are the error squared.
SE = np.sqrt(np.diag(covariance))
SE_A = SE[0]
SE_B = SE[1]

#Use the fit parameters to make a set of estamated y values 
#  from the fit equation. We can pass in the whole xdata array
#  and get out all the y value estimates at once using our 
#  function with our fit equation.  I called the new y-values
#  fit_y.
fit_y = RC_Charging(xdata, fit_A, fit_B)

#Plot the data (as dots) and the fit (as a line) to see if the equation
#  makes sense as a good fit.
plt.xlabel('Time in seconds')
plt.ylabel('Voltage (V)')
plt.plot(xdata, ydata, 'ro', label='data')
plt.plot(xdata, fit_y, 'b-', label='fit')
plt.legend()

#and print out our fit parameters and their uncertainties
print("  ")
print(F'The value of A is {fit_A:.5f} with standard error of {SE_A:.5f}.')
print(F'The value of B is {fit_B:.5f} with standard error of {SE_B:.5f}.')
print("  ")

#sometimes the curve fit routine throws a math warning, let the user know 
#  that the program ended and not to be upset about the warning
print("successful end of program, warning about overflow may follow")
print("  ")
\end{lstlisting}
\normalsize
It can do other languages, like the C language.

\small
\begin{lstlisting}[language=C] 
#include "ran0.h"
float ran0(int* idum)
{
	static float y,maxran,v[98];
	float dum;
	static int iff=0;
	int j;
	unsigned i,k;
	
	if((*idum < 0)||(iff == 0))
	{
		iff=1;
		i=2;
		maxran = RAND_MAX+1.0;
		*idum=1;
		for(j=1;j<97;j++)
		{
			dum=rand();
		}
		for(j=1;j<=97;j++)
		{
			v[j]=rand();
		}
		y=rand();
	}
	j=1+97.0*y/maxran;
	if((j>97)||(j<1))
	{
		cout << "This cannot happen." << endl;
		exit(0);
	}
	y = v[j];
	v[j] = rand();
	return y/maxran;
};
\end{lstlisting}
\normalsize

But if you don't need syntax highlighting (like if you will print in only black and white) you can also use a verbatim environment. In both cases, \LaTeX \enspace doesn't try to interpret your computer program code as \LaTeX \enspace commands. I want to make sure the code fits on the page, so for the verbatum envirnment example I am going to reduce the font size with a small environment as well. 
\begin{small}
\begingroup
\makeatletter
\@totalleftmargin=-1cm
	\begin{verbatim}
##############################################################
# Euler code to calculate a satellite orbit                  #
##############################################################
# Todd Lines
# 2022 04 19
##############################################################
# This code will calculate the path of an orbit given a 
#  specific velocity of the object and it's direction. 
#  It would probably be be better to have the whole calcuation 
#  go in reverse. But I didn't do that on purpose since we are
#  verifying that an elipse comes from a distance squared law
#  with this code.
##############################################################
# Load Libraries
import numpy as np                 # Numerical calculation library
import matplotlib.pyplot as plt    # Data plotting libary

##############################################################
# Initial conditions and physial setup constants
R=6371E3       # Radius of Earth in m
M=5.9E24       # Mass of plannet in Kg
m=0.020        # mass of our sattelite in kilograms
G=5.11E-11     # N*m**2/kg**2 Gravitational Constant
x0=2*R         # Initial x position in m (The Earth is at x = 0)
y0=0.5*R       # Initial y position in m (The Earth is at y = 0)
v0=6000.0       # Initial satellite velocity in m/s
thetadeg=100   # Satellite launch angle in degrees

##############################################################
## Set up the time steps and number of calcualtions
deltat=1.0        # Time steps of 0.01 seconds
ti=0              # starting at t=0
tf=300000         # final time in seconds
N=int((tf-ti)/deltat)  # calcualte how many time steps are in tf-ti seconds

##############################################################
# Preliminary calculations
pi=np.arccos(-1.0)         # calculate pi to machine percision
theta=thetadeg*pi/180      # calcualte theta in radians
vx1=v0*np.cos(theta)       # calculte the x component of the initial velocity
vy1=v0*np.sin(theta)       # calculte the y component of the initial velocity

##############################################################
# define and zero arrays
t=np.zeros((N))
x=np.zeros((N))
y=np.zeros((N))
vx=np.zeros((N))
vy=np.zeros((N))


##############################################################
## make an array of times to use
t=np.linspace(0,tf,num=N);
##############################################################
# Draw The Earth
r=R                               # distance from center in m
ThetaE0=0.00                      # initial angle in degrees
delta_ThetaE=5                    # change in angle in degrees
# calculate the number of points to use in our Earth drawing
Ne= int(360/delta_ThetaE)+1       # number of points
ye=np.zeros((Ne))
xe=np.zeros((Ne))
# Draw Circle to represent the Earth
ye[0]=0        # initial y positoin
xe[0]=r
thetaE=ThetaE0
for i in range (0,Ne):
thetaE=thetaE+delta_ThetaE
xe[i]=r*np.cos(thetaE*pi/180)
ye[i]=r*np.sin(thetaE*pi/180)
plt.plot(xe, ye)

##############################################################
# now perform an Euler's Method Calculation
# now recalling that vx(i) already has a cos(theta) in it,
# we can use this to calculate the x part of the resistive
# force and likewise use vy(i) in calculating the y part of
# the resistive force. No explicit calculation of theta is
# necessary this way, and we save lots of computation time.
x=np.zeros((N))
y=np.zeros((N))
vx=np.zeros((N))
vy=np.zeros((N))

# Put the starting point in the arrays
x[0]=x0                        # initial x position
y[0]=y0                        # initial y positoin
vx[0]=vx1                      # initial x velocity
vy[0]=vy1                      # initial y velocity, what we give it

# Now calculate all the other points of the satellite path
print('working')
for i in range (0,N-1):
r=np.sqrt(x[i]**2+y[i]**2)
if x[i] != 0:
theta=np.arctan(y[i]/x[i])
else:
theta=pi/2

if y[i]>0 and x[i]<0:
theta=pi+theta
if y[i]<0 and x[i]<=0:
theta=theta+pi
if y[i]<0 and x[i]>0:
theta=2*pi+theta
if r<R:
print("hit surface")
break   #hit surface

# Calculate the velocity and acceleration terms   
fx=vx[i]
gx=-np.cos(theta) * G*M/r**2
fy=vy[i]
gy=-np.sin(theta) * G*M/r**2

# Take an Euler step
x[i+1]=x[i]+deltat*fx
y[i+1]=y[i]+deltat*fy
vx[i+1]=vx[i]+deltat*gx
vy[i+1]=vy[i]+deltat*gy   
print('done')

# Now plot our satellite path 
plt.plot(x[:i],y[:i])
plt.show()

# And that is the end of the program
	
	\end{verbatim}
\endgroup
	
\end{small}