September 12, 2007
April 12, 2007
Circular Halo
I hardly noticed this without my Polaroid sunglasses on, as the colours were quite faint (and I’m fairly colour-blind), but the camera seemed to pick it out reasonably well; though its not as impressive as this one.
I tried to estimate the angular measurement of the halo and got my measurements wrong by quite a bit. From the sun to the start of the arc was the distance from my thumb to little finger tips with my hand held spread out at arms length. That’s 2 feet away from my eye and 8 inches across (I’ve measure that before and its easy to remember – all nice round numbers; no metric system centimeters here for estimating thank you – those units are too small.) This gives the halo an angular radius of about… (using a handy Common Lisp prompt)
CL-USER> (* (/ 180 pi) (atan 8/24))
18.434948654789338d0
that’s nowhere near the 22 degrees its supposed to be. But of course, holding your arm up at an angle reduces the distance between your eye and your hand that I’d measured before with my arm out horizontally. Re-measuring at around 45 degrees the hand-eye distance is down to 20″ and directly overhead 16″:
CL-USER> (* (/ 180 pi) (atan 8/20))
21.801408877810637d0
much better, phew.
April 5, 2007
Patterns
I think I got the matrix formula for this pattern from “Introduction to Dynamics” by Percival and Richards (its certainly the pattern I remember from the book cover)
You can get different patterns if you move the slider. Python gtk code…
import pygtk
pygtk.require('2.0')
import gtk, math
class Pattern(object):
x, y = 0, 0
points=[]
alpha = 76.11
def __init__(self, xsize, ysize):
self.xsize, self.ysize = xsize, ysize
window = gtk.Window(gtk.WINDOW_TOPLEVEL)
window.set_title("Pattern")
window.connect("destroy", lambda w: gtk.main_quit())
vbox = gtk.VBox(homogeneous=False, spacing=5)
window.add(vbox)
self.area = gtk.DrawingArea()
self.area.set_size_request(xsize, ysize)
self.area.connect("expose-event", self.area_expose_cb)
self.area.connect("size-allocate", self.calculate_points)
slider = gtk.Adjustment(value=self.alpha, lower=0.0,
upper=100.0, step_incr=0.01)
slider.connect("value-changed", self.slider_changed)
self.hscale = gtk.HScale(slider)
vbox.pack_start(self.area, expand=True, fill=True, padding=0)
vbox.pack_start(self.hscale, expand=False)
window.show_all()
def slider_changed(self, slider):
self.alpha = slider.get_value()
self.calculate_points()
def calculate_points(self, area=None, event=None):
x, y, self.xsize, self.ysize = self.area.get_allocation()
xscale = self.xsize/2.0
yscale = self.ysize/2.0
m = 52
angle = self.alpha*math.pi/180.0
c = math.cos(angle)
s = math.sin(angle)
self.points=[]
try:
for j in range(m):
x = 0
y = j/float(m)
for n in range(200):
w = x
x = x*c - (y - x*x)*s
y = w*s + (y - w*w)*c
if abs(x)>4 or abs(y)>4: raise StopIteration()
if x>1 or y>1: continue
self.points.append((int(x*xscale+xscale), int(y*yscale+yscale)))
except StopIteration:
self.area.queue_draw()
def area_expose_cb(self, area, event):
blue = self.area.get_colormap().alloc_color("#0000FF")
pointgc = self.area.window.new_gc()
pointgc.set_foreground(blue)
area.window.draw_point(pointgc, 20, 10)
self.area.window.draw_points(pointgc, self.points)
def main():
gtk.main()
return 0
if __name__ == '__main__':
Pattern(600, 500)
main()
November 16, 2006
Pattern Formation
cccc
cccc program from 'Models of Biological Pattern Formation' Hans Meinhardt
cccc
! gnuplot> set contour
! gnuplot> unset surface
! gnuplot> set style data lines
! gnuplot> set view 0,0
! gnuplot> set cntrparam levels 30
! gnuplot> splot 'ttt'
parameter (igridx = 31, igridy = 31)
real*8 g1(igridx,igridy), g2(igridx,igridy)
real*8 s1(igridx,igridy), s2(igridx,igridy)
real*8 vz(igridx,igridy)
real*8 dg1(30), dg2(30), dr(30), ds1(30), ds2(30)
kx = 30
ky = 30
qbb = 0.3
ginit = sqrt(1.d0 - qbb)
sinit = ginit
ta = 0.1
tb = 0.1
dia = 0.005
qd = 0.04
qe = 0.2
qa = 0.001
qaa = 0.001
cccc initialise the grid
do iy = 1, ky
do ix = 1, kx
g1(ix,iy) = ginit
g2(ix,iy) = ginit
s1(ix,iy) = sinit
s2(ix,iy) = sinit
c vz(ix,iy) = 1.d0 + rand()*0.0001
vz(ix,iy) = 1.d0 + rand()*0.0001
enddo
enddo
jy = (ky-1)/2 + 1
cccc initial perterbation
g1(kx,jy) = 1.1*ginit
do iter = 1, 5
dgc = 1.d0 - ta - 4.d0*dia
dsc = 1.d0 - qd - 4.d0*qe
drc = 1.d0 - tb
do it = 1, 500
do ix = 1, kx
g1(ix, ky+1) = g1(ix, ky) ! lower border
g2(ix, ky+1) = g2(ix, ky)
s1(ix, ky+1) = s1(ix, ky)
s2(ix, ky+1) = s2(ix, ky)
dg1(ix) = g1(ix, 1) ! upper border
dg2(ix) = g2(ix, 1)
ds1(ix) = s1(ix, 1)
ds2(ix) = s2(ix, 1)
enddo
do iy = 1, ky
g1(kx+1, iy) = g1(kx, iy) ! right border
g2(kx+1, iy) = g2(kx, iy)
s1(kx+1, iy) = s1(kx, iy)
s2(kx+1, iy) = s2(kx, iy)
enddo
do iy = 1, ky
g1l = g1(1,iy) ! left border
g2l = g2(1,iy)
s1l = s1(1,iy)
s2l = s2(1,iy)
cccc reaction and feedback
do ix = 1, kx
gf1 = g1(ix,iy)
gf2 = g2(ix,iy)
sf1 = s1(ix,iy)
sf2 = s2(ix,iy)
g1(ix,iy) = gf1*dgc+dia*
& (dg1(ix)+g1(ix,iy+1)+g1l+g1(ix+1,iy))+
& sf2*ta*vz(ix,iy)/(gf2**2+qbb)+qa
g2(ix,iy) = gf2*dgc+dia*
& (dg2(ix)+g2(ix,iy+1)+g2l+g2(ix+1,iy))+
& ta*sf1/(gf1**2+qbb)+qaa
s1(ix,iy) = sf1*dsc+qe*
& (ds1(ix)+s1(ix,iy+1)+s1l+s1(ix+1,iy))+qd*gf1
s2(ix,iy) = sf2*dsc+qe*
& (ds2(ix)+s2(ix,iy+1)+s2l+s2(ix+1,iy))+qd*gf2
g1l = gf1
g2l = gf2
dg1(ix) = gf1
dg2(ix) = gf2
s1l = sf1
s2l = sf2
ds1(ix) = sf1
ds2(ix) = sf2
enddo
enddo
enddo
enddo
cccc autocatalysis realised by inhibition of inhibition
C$$$ do iy = 1, igridy
C$$$ print '(31(g14.5,2x))', (s1(ix,iy), ix=1, igridx)
C$$$ enddo
do iy = 1, ky
do ix = 1, kx
print *, s1(ix,iy)
enddo
print *
enddo
end
November 7, 2006
Series Solution Method in maxima
Just a small script to illustrate series solution
/***************************************************************************
Solve the Hamiltonian
2 2 2
m q w p
H(p, q) = ------- + ---
2 2 m
by splitting to two coupled equations
d p
-- (q(t)) = -
dt m
d 2
-- (p(t)) = - m q w
dt
and using series method
****************************************************************************/
kill(p,q,t,m,w)$
eq1: diff(q(t), t)= p(t)/m;
eq2: diff(p(t), t)= -m*w^2*q(t);
assume(w>0)$
atvalue(p(t), t=0, 0)$
atvalue(q(t), t=0, 1)$
ans: desolve([eq1, eq2], [q(t), p(t)]);
/* And now do it just using power series */
kill(a,b,p,q,t,w,m)$
p:taylor(sum(a[i]*t^i,i,0,6),t,0,6)$
q:taylor(sum(b[i]*t^i,i,0,6),t,0,6)$
/* Note rearrange eqn=0 */
eq1:diff(q,t)-p/m$
eq2:diff(p,t)+m*w^2*q$
coeffs1:makelist(ratcoef(eq1,t,i),i,0,6)$
coeffs2:makelist(ratcoef(eq2,t,i),i,0,6)$
eqns:append(cons(subst(0,t,p) = 0,coeffs1),cons(subst(0,t,q) = 1,coeffs2))$
ab:append(makelist(a[i],i,0,6),makelist(b[i],i,0,6))$
ans:solve(eqns,ab)$
trunc(subst(ans[1],p));
trunc(subst(ans[1],q));
or inside an emaxima session:
% -*-EMaxima-*-
\documentclass{article}
\usepackage{a4wide}
\usepackage{emaxima}
\begin{document}
\section*{Series Solution of PDE}
$$
H(p,q) = {p^2\over 2 m} + {m\over 2} w^2 q^2
$$
$$
\frac{dq}{dt}=\frac{\partial H}{\partial p} = \frac{p}{m} \qquad{\rm and}\qquad
\frac{dp}{dt}=-\frac{\partial H}{\partial q} = -m\omega^2q
$$
\begin{maximasession}
kill(p,q,t,m,w)$
eq1: diff(q(t), t)= p(t)/m;
eq2: diff(p(t), t)= -m*w^2*q(t);
assume(w>0)$
atvalue(p(t), t=0, 0)$
atvalue(q(t), t=0, 1)$
ans: desolve([eq1, eq2], [q(t), p(t)]);
\maximaoutput*
\end{maximasession}
To do that with series expansions only:
\begin{maximasession}
p: taylor(sum(a[i]*t^i,i,0,inf),t,0,6)$
q: taylor(sum(b[i]*t^i,i,0,inf),t,0,6)$
eq1: diff(q, t) - p/m$
eq2: diff(p, t) + m*w^2*q$
coeffs1: makelist (ratcoeff(eq1, t, i), i , 0 , 6)$
coeffs2: makelist (ratcoeff(eq2, t, i), i , 0 , 6)$
/* add in boundary conditions p=0 and q=1 at t=0 */
eqns: append(cons(subst(0,t,p)=0, coeffs1), cons(subst(0,t,q)=1, coeffs2))$
coeffs:append(makelist(a[i],i,0,6), makelist(b[i],i,0,6))$
ans: solve(eqns, coeffs)$
trunc(subst(ans[1], p));
trunc(subst(ans[1], q));
\maximaoutput*
\end{maximasession}
\end{document}
November 4, 2006
emaxima
Usually I run maxima just from the command line as rmaxima but occasionally it is useful to run it inside emacs in combination with Tex.
My .emacs has
;;;
;;; Maxima
;;;
(add-to-list 'load-path "/home/r/src/lisp/maxima/interfaces/emacs/emaxima")
(autoload 'maxima "maxima" "Running Maxima interactively" t)
(autoload 'maxima-mode "maxima" "Maxima editing mode" t)
(autoload 'emaxima-mode "emaxima" "eMaxima" t)
(add-hook 'inferior-maxima-mode-hook
(lambda () (local-set-key [tab] 'inferior-maxima-complete)))
(add-hook 'inferior-maxima-mode-hook
(lambda () (set (make-local-variable 'comint-prompt-read-only) t)))
and an example of a tex file that includes some maxima workings…
% -*-EMaxima-*-
\documentclass{article}
\usepackage{a4wide} % bit more room!
\usepackage{amsmath} % used to get e.g. \boxed{...}
\usepackage[breqn]{emaxima}
\begin{document}
\section*{Laplace Transform for Solving Integral Equation}
% C-c + to go to next cell.
% C-u C-c C-u a updates all the cells without asking - maxima output
% C-u C-c C-u A updates all the cells without asking - TeX output
Solve $$\int_0^t cosh(a x) f(t-x)\,dx + b f(t)=t^3$$
\begin{maximasession}
'integrate(cosh(a*x)*f(t-x), x, 0, t) + b*f(t)=t**3;
laplace(%,t,s);
res: linsolve([%],['laplace(f(t),t,s)]);
\end{maximasession}
Use {\tt dpart} to put a box around which part you will get
when using {\tt part}
% use C-c C-u C to get the tex output for just this one
% C-c C-d to delete the output
\begin{maxima}[laplace1]
dpart(res[1], 2);
ilt(part(res[1], 2),s,t);
\end{maxima}
\end{document}
and a slightly more complex example…
% -*-EMaxima-*-
\documentclass{article}
%\usepackage{a4wide} % bit more room!
\usepackage{amsmath} % used to get e.g. \boxed{...}
\usepackage[breqn]{emaxima}
\usepackage{vpage}
\setpapersize{Afour}
\setmarginsrb{1in}{1in}{1in}{1in}{0pt}{0mm}{0pt}{0mm}
\begin{document}
\section*{Perturbed DE Solution}
Solve $u''+u=\epsilon u^2$ with $u(0)=1/r$ and $u'(0) = 0$ by looking
for solutions of the form
$$
u(t)=u_0(t)+\epsilon u_1(t)+\epsilon^2u_2(t)+\ldots
$$
and initial conditions $u_0 = 1/r$, $u'_0=0$,
$u_1=0$, $u'_1=0$, $u_2=0$, $u'_2=0$ at $t=0$
\begin{maximasession}
eqn_init: 'diff(u,t,2)+u=epsilon*u^2;
v: u0 + epsilon*u1 + epsilon^2*u2;
depends([u0,u1,u2],t);
eq: subst(v,u,eqn_init);
eq: ev(eq, diff);
eq: ratexpand(eq);
equ0: ev(eq, epsilon=0);
solveu0: ode2(equ0, u0, t);
solveu0: ic2(solveu0, t=0, u0=1/r, 'diff(u0,t)=0);
uu0: rhs(%)$
eq1: coeff(eq, epsilon, 1);
eq1: ev(eq1, u0=uu0);
solveu1: ode2(eq1, u1, t);
solveu1: ic2(solveu1, t=0, u1=0, 'diff(u1,t)=0);
uu1: rhs(%)$
eq2: coeff(eq, epsilon, 2);
eq2: ev(eq2, u0=uu0, u1=uu1);
eq2: ratexpand(eq2);
solveu2: ode2(eq2, u2, t);
solveu2: ic2(solveu2, t=0, u2=0, 'diff(u2,t)=0);
uu2: rhs(%)$
u_final: ev(v, u0=uu0, u1=uu1, u2=uu2);
trunc(%);
\end{maximasession}
\end{document}
February 23, 2005
Snowbow
ok, I made that name up. It seem there is officially no such thing as a snowbow, but I think its a good name (as it was starting to snow when I took this picture).
So why is it different from a rainbow? Well, when you see a rainbow the sun is always behind you and the sunlight is being refracted from water droplets back towards you – but in the picture the sun is infront of me and I’m leaning out of my window and looking up (you can just see the guttering at the eaves of the house at the top of the picture).
More correctly its called a ‘Circumzenithal Arc’ according to this excellent rainbow site.
June 8, 2004
Transit of Venus
In a very peculiar coincidence, a few months ago I was talking to a neighbour about a little lane that runs off from the road a hundred yards or so up the hill from where I live, and they mentioned there was a blue plaque to commemorate an astronomer at the bottom. Intrigued, I went and had a look and sure enough there it was on a garden wall.
Looking up one the net for information I could not believe that there was shortly going to be another transit, the first for over a hundred years. So I had to take the day off work and try and record the event in honour of William.
When I was a teenager I had borrowed a small telescope and was lucky enough to have the right positioning to project the sun onto the back white wall of my bedroom. With thick curtains darkening the room and pegged together around the scope tube I achieved an image of the sun several feet across and, apart from some shimmering as the hot air rose up the walls of the house, this gave some great views of sunspots.
The weather was blazing sunshine and with my Meade ETX, some cardboard a sheet of white paper and a certain amount of swearing my contraption held steady enough for a few pictures.
September 16, 2001
Countdown numbers game
No this isn’t physics, but I suppose it charitably comes under the heading of numerical methods ;) …
I wrote the following (probably badly written) program after watching Carol Vorderman struggle to find solutions for the numbers game on the Countdown TV program.
I can’t remember why I wrote it in Fortran 95 ?!? and I didn’t test it much, just on a couple of sets of numbers I don’t think she (or I!) could solve by hand in the 30 secs the contestants are given.
The problem is – given a set of 6 random integers try to arrive at a randomly chosen number just using +,-,*,/ (although, thinking about it, I’m not sure integer division is allowed)
For the second test case (try to use 50, 9, 9, 4, 7 and 10 to make 971) the program printed
***SOLUTION***
9 * 4 = 36
7 * 10 = 70
50 - 36 = 14
70 * 14 = 980
980 - 9 = 971
***SOLUTION***
9 * 7 = 63
50 + 63 = 113
113 - 4 = 109
9 * 109 = 981
981 - 10 = 971
at least 2 solutions
and the code COUNTDOWN.F90 is…
module solutions
integer, parameter::ncards=6
integer::nsols
integer::sol(1000,ncards*2)
end module solutions
program countdown
use solutions
integer::numbers(ncards)
integer::solunum(ncards)
integer::solunum1(ncards)
integer::ntarget
character::soluop(ncards)
numbers=(/50, 7, 5, 1, 8, 3/)
ntarget=684
numbers=(/50, 9, 9, 4, 7, 10/)
ntarget=971
solunum=0
solunum1=0
soluop=' '
nsols=1
call pair(numbers, ncards, solunum, solunum1, soluop, ncards, ntarget)
print *, 'at least ',nsols-1, ' solutions'
end program countdown
integer function combine(n, m, op)
character::op
integer::n, m
select case (op)
case ('*')
combine=n*m
case ('+')
combine=n+m
case ('-')
combine=n-m
case ('/')
combine=n/m
if (mod (n,m) /= 0) combine=0
case default
combine=0
end select
if (combinesolsize*2) then
print *, '***SOLUTION***'
do m=solsize, 2, -1
n1=solunum(m)
n2=solunum1(m)
select case (soluop(m))
case ('*')
n3=n1*n2
case ('+')
n3=n1+n2
case ('-')
n3=n1-n2
case ('/')
n3=n1/n2
end select
print '(i4, 2x, a, i4, 2x, a, i4)', n1,soluop(m), n2, '=', n3
enddo
sol(nsols,:)=nums
nsols=nsols+1
return
endif
enddo
end subroutine printsolution
October 27, 1995
Teaching MSc Numerical Methods
Around this time I taught a few groups of MSc students some basic numerical methods – mainly to support another class in soliton theory.
One of the first programs was to plot optical reflectance against angle at a metal glass interface – this is interesting because at certain angles you will get a surface wave traveling in the metal film.
Mostly though this was just to get students comfortable using 2D plotting NAG routines…
(more…)
November 28, 1994
Graphics and Plotting during my Ph.D
It seems like a considerable chunk of my time during my PhD was devoted to getting half decent graphics output. When I started everyone was using NAG which was fast but very basic. I got handed a three folders full of GINO manuals and asked to write something that would plot in colour on the new HP Deskjet printers.
I remember being very fed up that the department had paid several hundred pounds for a graphics library that, although it contained routines for fitting curves to datapoints, wouldn’t scale the graph axes to cope with that – you had to do that yourself!
CCCC ===================================================================
CCCC Converts the output from GINO type plot files for GNUPLOT to use
CCCC ===================================================================
CCCC -------------------------------------------------------------------
CCCC Usage GCON "infile" "outfile"
CCCC Version 1.0
CCCC Created 28/11/94[1.0]
CCCC -------------------------------------------------------------------
PROGRAM Convert
PARAMETER (mdim = 150)
DOUBLE COMPLEX pcmplx(mdim, mdim)
REAL pcoor(mdim, mdim)
CHARACTER*50 buf
CHARACTER*50 title
CHARACTER*1 choice
INTEGER*4 count
INTEGER*2 stat, i
count = NARGS()
IF (count .NE. 3) THEN
PRINT *, 'SYNTAX: GCON infile outfile '
STOP
END IF
CALL GETARG(1, buf, stat)
OPEN (10, file = buf(1:stat), STATUS='OLD')
CALL GETARG(2, buf, stat)
OPEN (11, file = buf(1:stat))
READ(10,'(A)') title
READ(10,*) nx,xmin,xmax,ymin,yinc
do j=1,mdim
read(10,*,end=10) (pcmplx(i,j),i=1,nx)
enddo
print *,'****array out of bounds, increase mdim****'
stop
10 CONTINUE
ny = j - 1
PRINT *,' Real (R), Imaginary (I), Absolute Squared (A) ?'
READ '(A)', choice
IF (choice.eq.'R' .OR. choice.eq.'r') THEN
DO i = 1, nx
DO j = 1, ny
pcoor(i,j) = REAL(pcmplx(i,j))
END DO
END DO
ELSEIF (choice.eq.'I' .OR. choice.eq.'i') THEN
DO i = 1, nx
DO j = 1, ny
pcoor(i,j) = AIMAG(pcmplx(i,j))
END DO
END DO
ELSEIF (choice.eq.'A' .OR. choice.eq.'a') THEN
DO i = 1, nx
DO j = 1, ny
pcoor(i,j) = CABS(pcmplx(i,j))**2
END DO
END DO
ENDIF
xinc = (xmax - xmin) / nx
DO j = 1, ny
DO i = 1, nx
WRITE(11,20) xmin + xinc*(i-1), ymin+(j-1)*yinc, pcoor(i,j)
cccc WRITE(6,20) xmin + xinc*(i-1), ymin+(j-1)*yinc, pcoor(i,j)
END DO
WRITE (11,'(A)')
END DO
20 FORMAT (1x, 3(E12.6,1x))
CLOSE(10)
CLOSE(11)
END
April 30, 1993
BSc Fortran Numerics Course – Phonon Dispersion Curves
Amongst the programs were those to calculate phonon dispersion curves and had a note…
“The frequency of oscillation of an atom in a crystal can be calculated from an equation involving a dynamical matrix which incorportes the mass of the atom and the forces on it.
The frequency of oscillation in three crystal directions [100], [110] and [111] are calculated for iron. The program includes data such as the atomic weights and elastic bulk modulus and requests input of the number of atoms in the unit cell which is 2 for BCC and 4 for FCC.”
CCCC*******************************************************************
CCCC Calculation of phonon frequencies travelling through a
CCCC iron crystal lattice, for [100] [110] and [111] directions
CCCC*******************************************************************
C NOTES: User input determines number of atoms per unit cell
C BCC=2 FCC=4
C
C VARIABLES:
C RHO - Density
C K - Bulk modulus
C N - Number of atoms per unit cell
C KX,KY,KZ - Wave direction components
C A - Unit cell size
C PIOA - Pi over A
C STEP - Stepsize for increment to Kx
C AT - Atomic weight
C GAMMAM - Force constant divided by mass
C
C SUBROUTINES:
C F02ABF - NAG Routine calculates the eigenvalues and eigenvectors
C of a symmetric matrix.
C PH - Programs the dynamical matrix and calls F02ABF.
C
CCCC*******************************************************************
REAL*8 RHO,K,KX,KY,KZ,A,PIOA,STEP,AT,GAMMAM,R
CCCC Data for Iron
DATA AT,RHO,K,R/55.847,7860.0D0,170.0D9,0.8D0/
CCCC Read number of atoms in the unit cell
PRINT *,'ENTER THE NUMBER OF ATOMS IN UNIT CELL '
READ *, N
CCCC Calculate unit cell size
A = (N*AT/(1000.0*6.023D23*RHO))**(1.0/3.0)
GAMMAM = K/(3.0D0*(4.0D0+4.0D0*R)*RHO*A*A)
PIOA = 4.0d0*DATAN(1.0d0)/A
STEP = PIOA/12.0D0
PRINT *,' KX [100] FREQUENCIES'
DO 10 KX = 0.0D0, PIOA, STEP
CCCC Calculated wave frequencies for [100] propagation direction
KY = 0.0D0
KZ = 0.0D0
CALL PH(KX,KY,KZ,AT,GAMMAM,R,A)
10 CONTINUE
PRINT *,' KX [110] FREQUENCIES'
DO 20 KX = 0.0D0, PIOA, STEP
CCCC Calculated wave frequencies for [110] propagation direction
KY = KX
KZ = 0.0D0
CALL PH(KX,KY,KZ,AT,GAMMAM,R,A)
20 CONTINUE
PRINT *,' KX [111] FREQUENCIES'
DO 30 KX = 0.0D0, PIOA, STEP
CCCC Calculated wave frequencies for [111] propagation direction
KY = KX
KZ = KX
CALL PH(KX,KY,KZ,AT,GAMMAM,R,A)
30 CONTINUE
END
SUBROUTINE PH(KX,KY,KZ,AT,GAMMAM,R,A)
CCCC*******************************************************************
C Finds the eigenvalues - the wave frequencies for a specified
C wave direction
CCCC*******************************************************************
C VARIABLES:
C A - Unit cell size
C AT - Atomic Weight
C GAMMAM - Force constant divided by mass
C R - Force constant ratio
C COEFF - Coefficient
C DELTA - Constant
C D(,) - Dynamical Matrix
C KX,KY,KZ - Wave direction components
C
C ADDITIONAL VARIABLES FOR NAG ROUTINE:
C IA - Specifies the first dimension of D(,)
C ORDER - Specifies the order of the matrix
C IFAIL - ERROR INDICATOR
C EIGVAL() - Contains eigenvalues
C EIGVEC(,) - Contains eigenvectors
C IV - First dimension of eigenvector array
C E - Work space
C
C SUBROUTINES:
C F02ABF - NAG Routine calculates the eigenvalues and eigenvectors
C of a symmetric matrix.
CCCC*******************************************************************
EXTERNAL F02ABF
REAL*8 KX,KY,KZ,A,GAMMAM,D(3,3),EIGVAL(3),COEFF,
* EIGVEC(3,3),DELTA,R,AT,E(3)
INTEGER ORDER
CCCC Set values for NAG routine
IA = 3
ORDER = 3
IFAIL = 0
IV = 3
CCCC Calculate Constants
COEFF = 8.0D0/3.0D0
DELTA = COEFF+2.0D0*R - COEFF*DCOS(KX*A)*DCOS(KY*A)*DCOS(KZ*A)
CCCC Calculated dynamical matrix elements
D(1,1) = DELTA-2.0D0*R*DCOS(2.0*KX*A)
D(2,1) = COEFF*DSIN(KX*A)*DSIN(KY*A)*DCOS(KZ*A)
D(3,1) = COEFF*DSIN(KX*A)*DCOS(KY*A)*DSIN(KZ*A)
D(1,2) = D(2,1)
D(2,2) = DELTA-2.0D0*R*DCOS(2.0*KY*A)
D(3,2) = COEFF*DCOS(KX*A)*DSIN(KY*A)*DSIN(KZ*A)
D(1,3) = D(3,1)
D(2,3) = D(3,2)
D(3,3) = DELTA-2.0D0*R*DCOS(2.0*KZ*A)
CALL F02ABF(D,IA,ORDER,EIGVAL,EIGVEC,IV,E,IFAIL)
IF (IFAIL.NE.0) THEN
PRINT *,'ERROR IN NAG ROUTINE IFAIL = ',IFAIL
STOP
ENDIF
CCCC Multiply eigenvalues by force constant divided by mass
PRINT '(1X,4(D12.5,1X))',KX,(EIGVAL(I)*GAMMAM,I=1,3)
END
I think the programs were part of some coursework I did during my BSc in Computational Physics at Salford University.
It must have been after my stint at BNFL because I was using raw TeX as I thought Latex looked too babyish. I’d put together some basic macros…
\\font\\ninerm=cmr9
\\font\\twelverm=cmr12
\\font\\seventeenrm=cmr17
\\font\\smallcaps=cmcsc10
\\font\\bigtitle=cmss12
\\font\\part=cmbx12
\\parindent 0mm
\\parskip 1em
\\baselineskip=1.2em
\\def\\E{{\\bf E}}
\\def\\B{{\\bf B}}
\\def\\H{{\\bf H}}
\\def\\P{{\\bf P}}
\\def\\Pl{{\\bf P}_L}
\\def\\Pnl{{\\bf P}_{NL}}
\\def\\D{{\\bf D}}
\\def\\r{{\\bf r}}
\\def\\del{\\nabla}
\\def\\intf{\\intop^{+\\infty}_{-\\infty}}
\\def\\chione{\\chi^{(1)}}
\\def\\chithree{\\chi^{(3)}}
\\def\\muo{\\mu_0}
\\def\\ko{k_0}
\\def\\wo{w_0}
\\def\\betao{\\beta_0}
\\def\\epsilono{\\epsilon_0}
\\def\\Eh{\\tilde E}
\\def\\Ah{\\tilde A}
\\def\\Eb{{\\overline E}}
\\def\\nb{{\\overline n}}
\\def\\betab{{\\overline\\beta}}
\\def\\pone#1#2{{\\partial#1\\over\\partial#2}}
\\def\\ptwo#1#2{{\\partial^2 #1\\over\\partial #2^2}}
\\def\\Tos{T_0^2}
\\def\\sqp{\\sqrt{P_0}}
\\def\\fig#1{{\\bf Figure #1. }}
%\\def\\section#1{{\\vskip 5mm\\seventeenrm #1\\vskip 1em}}
\\def\\section#1{\\vfill\\eject{\\seventeenrm #1\\vskip 1em}}
\\def\\subsection#1{\\vskip 5mm{\\bigtitle #1\\vskip 1em}}
\\def\\ve{\\vfill\\eject}
\\def\\fortran{{\\tt FORTRAN }}
\\def\\pagenumbers{\\footline={\\hss\\tenrm\\folio\\hss}}
\\def\\boxit#1{\\vbox{\\hrule\\hbox{\\vrule\\kern3pt
\\vbox{\\kern3pt#1\\kern3pt}\\kern3pt\\vrule}\\hrule}}
\\def\\uncatcodespecials{\\def\\do##1{\\catcode`##1=12 }\\dospecials}
\\newcount\\lineno % the number of file lines
{\\catcode`\\`=\\active \\gdef`{\\relax\\lq}}
\\def\\myupverbatim{\\tt \\lineno=0
\\def\\par{\\leavevmode\\endgraf} \\catcode`\\`=\\active
\\obeylines \\uncatcodespecials \\obeyspaces
\\everypar{\\advance\\lineno by1 \\llap{\\fiverm\\ \\ }}}
\\def\\setupverbatim{\\tt \\lineno=0
\\def\\par{\\leavevmode\\endgraf} \\catcode`\\`=\\active
\\obeylines \\uncatcodespecials \\obeyspaces
\\everypar{\\advance\\lineno by1 \\llap{\\fiverm\\the\\lineno\\ \\ }}}
{\\obeyspaces\\global\\let =\\ } % let active space = control space
\\def\\listing#1{\\par\\parskip -0.3em\\begingroup\\setupverbatim\\input#1 \\endgroup
\\parskip 1em}
\\def\\results#1{\\par\\parskip -0.3em\\begingroup\\myupverbatim\\input#1 \\endgroup
\\parskip 1em}
and I was using the fabulous GLE (Graphics Language Editor) by Chris Pugmire for my graphs (sigh, I still dont think anything looks as good, it had a cool trick to use TeX to do the axes captions). A typical gle input file was
size 21 29
amove 3 14
begin scale 0.6 0.6
begin graph
nobox
size 15 15
title "Current in a Diode Circuit" font texcmr
xtitle "Voltage /V" font texcmtt
ytitle "Current /Amps" font texcmtt
data vidata.dat
d1 lstyle 1 smooth
end graph
rmove 0 -14
begin graph
nobox
size 15 15
title "Solution Times for NAG and Bisection" font texcmr
xtitle "Voltage /V" font texcmtt
ytitle "Time for Solution /Secs" font texcmtt
data tdata.dat
yaxis max 1e-2
d1 lstyle 1 smooth key NAG
d2 lstyle 2 smooth key Bisection
end graph
end scale
April 4, 1993
Final Year BSc Project
I did some initial work on the effects of soliton propagation in birefringent fibres around this time writing a program I later cleaned up a bit for my PhD.
There were a few dozen (literally) versions of this lying around on my HDD, back then you had to write a new program and compile it for each new set of parameters you were investigating – which meant leaving lots of versions of code around.
Its one of the things that drove me to buy Mathematica 2.0 (I think it was). I developed a system for splicing the output from mathematica into a fortran program and running that – but it wasn’t general enough to be useful. The mathematica syntax, which I struggled using, led me to Scheme and finally to Common Lisp – I wish I could have used it back then!
Essentially I wrote a program to solve:
![{\partial u\over\partial\xi}=-\delta{\partial u\over\partial\tau}+{i\over 2}{\partial^2 u\over\partial\tau^2}+i\left[|u|^2+{2\over 3}|v|^2\right] u](http://s1.wordpress.com/latex.php?latex=%7B%5Cpartial+u%5Cover%5Cpartial%5Cxi%7D%3D-%5Cdelta%7B%5Cpartial+u%5Cover%5Cpartial%5Ctau%7D%2B%7Bi%5Cover+2%7D%7B%5Cpartial%5E2+u%5Cover%5Cpartial%5Ctau%5E2%7D%2Bi%5Cleft%5B%7Cu%7C%5E2%2B%7B2%5Cover+3%7D%7Cv%7C%5E2%5Cright%5D+u&bg=ffffff&fg=000000&s=3 "{\partial u\over\partial\xi}=-\delta{\partial u\over\partial\tau}+{i\over 2}{\partial^2 u\over\partial\tau^2}+i\left[|u|^2+{2\over 3}|v|^2\right] u")
![{\partial v\over\partial\xi}=+\delta{\partial v\over\partial\tau}+{i\over 2}{\partial^2 v\over\partial\tau^2}+i\left[|v|^2+{2\over 3}|u|^2\right] v](http://s2.wordpress.com/latex.php?latex=%7B%5Cpartial+v%5Cover%5Cpartial%5Cxi%7D%3D%2B%5Cdelta%7B%5Cpartial+v%5Cover%5Cpartial%5Ctau%7D%2B%7Bi%5Cover+2%7D%7B%5Cpartial%5E2+v%5Cover%5Cpartial%5Ctau%5E2%7D%2Bi%5Cleft%5B%7Cv%7C%5E2%2B%7B2%5Cover+3%7D%7Cu%7C%5E2%5Cright%5D+v&bg=ffffff&fg=000000&s=3 "{\partial v\over\partial\xi}=+\delta{\partial v\over\partial\tau}+{i\over 2}{\partial^2 v\over\partial\tau^2}+i\left[|v|^2+{2\over 3}|u|^2\right] v")
where u and v are different modes of the fibre. Fortran code…
(more…)
June 19, 1992
Fractals
This is a historical blog entry of some old Fortran files I found on my HDD…
I think this was my first attempt to draw fractals, unfortunately now (in 2007) I no longer have the compiler to try reproducing the graphs.
PROGRAM FRACTAL
INTEGER NTRANS,XSCALE,YSCALE,XOFFSET,YOFFSET,IY,IX
REAL A,B,C,D,E,F,P,X,Y,NEWX,NEWY
DOUBLE PRECISION PK,N
CHARACTER*9 NAME
DIMENSION A(100),B(100),C(100),D(100),E(100),F(100),
* P(100)
PRINT *,'ENTER FILENAME FOR DATA 5CHAR.DAT'
READ 100,NAME
OPEN (5,FILE=NAME,STATUS='OLD')
CALL ZEROISE(A,B,C,D,E,F,P)
CALL GETDATA(A,B,C,D,E,F,P,NTRANS,XSCALE,YSCALE)
CALL SETSCREEN(X,Y,XOFFSET,YOFFSET)
CALL DATE_TIME_SEED@
1 DO 10 N=1,200000
PK=RANDOM()
DO 5 M=NTRANS,1,-1
IF (PK.LE.P(M)) K=M
5 CONTINUE
NEWX=A(K)*X+B(K)*Y+E(K)
NEWY=C(K)*X+D(K)*Y+F(K)
X=NEWX
Y=NEWY
IX=INT(X*XSCALE+XOFFSET)
IY=INT(Y*YSCALE+YOFFSET)
CALL GET_KEY1@(J)
IF (J.NE.0) GOTO 12
IF (N.LT.10) GOTO 10
CALL SET_PIXEL@(IX,IY,14)
10 CONTINUE
12 IF (J.EQ.73) THEN
XSCALE=XSCALE+10
YSCALE=YSCALE+10
ENDIF
IF (J.EQ.79) THEN
XSCALE=XSCALE-10
YSCALE=YSCALE-10
ENDIF
IF (J.EQ.328) YOFFSET=YOFFSET-10
IF (J.EQ.336) YOFFSET=YOFFSET+10
IF (J.EQ.331) XOFFSET=XOFFSET-10
IF (J.EQ.333) XOFFSET=XOFFSET+10
CALL CLEAR_SCREEN@
IF (J.NE.81) GOTO1
CALL TEXT_MODE@
100 FORMAT (A)
END
SUBROUTINE ZEROISE(A,B,C,D,E,F,P)
REAL A,B,C,D,E,F,P,PK
DIMENSION A(100),B(100),C(100),D(100)
* ,E(100),F(100),P(100)
DO 1000 I=1,100
A(I)=0.0
B(I)=0.0
C(I)=0.0
D(I)=0.0
E(I)=0.0
F(I)=0.0
P(I)=0.0
1000 CONTINUE
END
SUBROUTINE GETDATA(A,B,C,D,E,F,P,NTRANS,XSCALE,YSCALE)
INTEGER NTRANS,XSCALE,YSCALE
REAL A,B,C,D,E,F,P,PK,PT
DIMENSION A(100),B(100),C(100),D(100)
* ,E(100),F(100),P(100)
READ (5,*) NTRANS
READ (5,*) XSCALE,YSCALE
PT=0
DO 2000 I=1,NTRANS
READ (5,*) A(I),B(I),C(I),D(I),E(I),F(I),PK
PT=PT+PK
P(I)=PT
2000 CONTINUE
END
SUBROUTINE SETSCREEN(X,Y,XOFFSET,YOFFSET)
INTEGER XOFFSET,YOFFSET
REAL X,Y
X=0.0
Y=0.0
XOFFSET=220
YOFFSET=0
CALL VGA@
CALL CLEAR_SCREEN@
END
and there were several datafiles of parameters which I think drew fern patterns etc.
4
50 40
0,0,0,0.16,0,0,0.01
0.2,-0.26,0.23,0.22,0,1.6,0.07
-0.15,0.28,0.26,0.24,0,0.44,0.07
0.85,0.04,-0.04,0.85,0,1.6,0.85
4
800 800
0 0 0 0.5 0 0 0.05
0.1 0 0 0.1 0 0.2 0.15
0.42 -0.42 0.42 0.42 0 0.2 0.4
0.42 0.42 -0.42 0.42 0 0.2 0.4
3
500 500
0.5 0 0 0.5 0 0 0.33
0.5 0 0 0.5 1 0 0.33
0.5 0 0 0.5 0.5 0.5 0.34
There was also a very primitive mandelbrot program which I managed to make work under wine and take a screenshot!
PROGRAM MANDEL
DIMENSION ZRL(60),ZPX(60)
SCALEX=266.666
SCALEY=200.0
LOFFX=320
LOFFY=240
CALL VGA@
CALL CLEAR_SCREEN@
DO 30 CPX=1.2,-1.2,-5.0E-3
DO 20 CRL=-1.2,1.2,3.75E-3
ZRL(1)=0.0
ZPX(1)=0.0
ZRL(2)=ZRL(1)*ZRL(1)-ZPX(1)*ZPX(1)+CRL
ZPX(2)=2.0*ZRL(1)*ZPX(1)+CPX
IF (ABS(ZRL(2)).GT.1E10) GOTO 11
IF (ABS(ZPX(2)).GT.1E10) GOTO 11
ZRL(3)=ZRL(2)*ZRL(2)-ZPX(2)*ZPX(2)+CRL
ZPX(3)=2.0*ZRL(2)*ZPX(2)+CPX
IF (ABS(ZRL(3)).GT.1E10) GOTO 11
IF (ABS(ZPX(3)).GT.1E10) GOTO 11
DO 10 I=3,50,1
ZRL(I+1)=ZRL(I)*ZRL(I)-ZPX(I)*ZPX(I)+CRL
ZPX(I+1)=2.0*ZRL(I)*ZPX(I)+CPX
DET=((ZRL(I)-ZRL(I-1))-(ZRL(I-1)-ZRL(I-2)))/2.0
DET1=((ZPX(I)-ZPX(I-1))-(ZPX(I-1)-ZPX(I-2)))/2.0
C PRINT *,I,DET,DET1
IF (ABS(DET).GT.1.0 .OR. ABS(DET1).GT.1.0) THEN
CALL SET_PIXEL@(INT(LOFFX+CRL*SCALEX)
* ,INT(LOFFY-CPX*SCALEY),15)
GOTO 20
ENDIF
10 CONTINUE
11 CALL SET_PIXEL@(INT(LOFFX+CRL*SCALEX)
* ,INT(LOFFY-CPX*SCALEY),0)
20 CONTINUE
30 CONTINUE
END
February 15, 1991
BSc Fortran Numerics Course
This a retro blog entry for the first ever fortran programming courses at Salford (taught by Dr Keeler iirc). I’ve had a FORTRAN directory kicking around on my machine for years containing all my old files so I thought I’d clean things up and create a blog entry for posterity.
How do I know the date so accurately?… well some of the files were log files generated from the compiles like…
Prospero Fortran Compiler - iid 2.1 Date: 1991 February 15 16:41:17.28
Compilation of: A:\BESSEL\DRUMDRW.FOR
Options: A1 B0 C0 E3 G1 H0 I1 L0 M0 N3 O1 Q0 R0 S1 T0 U1 V1 Y1
5 === 5 === 0 INCLUDE C:\PUBLIC\PROFOR\INCLUDE\F77PC.FOR
Ah yes, Prospero fortran, and we wrote our code on 3 1/2 inch floppies which we had to carry around with us. I think we used Prospero, even though the Salford
compiler FTN77 was available because the graphics were much easier to produce as they were built in. I remember it being a shock when I went over to FTN77 which found lots of bugs in my code that Prospero had allowed, guess that is why I went to work for Salford all those years later.
We used Fortran as just a calculator really e.g.
CCCC ******************************************************************
CCCC ** PRINTS THE VALUE OF THE BESSEL FUNCTION FOR X=0-20 **
CCCC ** AND ALLOWS BOUNDARIES FOR SOLUTION **
CCCC ** BY BISECTION TO BE FOUND **
CCCC ******************************************************************
REAL X,Y,TOL,BESS
PRINT99,'X','Y'
DO 10 X=0,20
Y=BESS(X,1.0E-5)
CALL DISP (X,Y)
10 CONTINUE
99 FORMAT (2(10X,A))
END
REAL FUNCTION BESS(X,TOL)
REAL X,TOL,OLD,TERM,J
INTEGER N
TERM=1
J=1
N=0
DO WHILE (ABS(TERM-OLD).GE.TOL)
OLD=TERM
N=N+1
TERM=-(OLD*X*X)/(4.0*N*N)
J=J+TERM
END DO
120 BESS=J
END
SUBROUTINE DISP(X,Y)
REAL X,Y
INTEGER XN,YN,OLDX,OLDY
IF (X.EQ.0) THEN
OLDX=20
OLDY=175
ENDIF
XN=INT(X*30+20)
YN=INT(175-Y*170)
PRINT*,X,Y
OLDX=XN
OLDY=YN
END
or solving simple equations that we’d use some symbolic maths package these days to solve…
CCCC ******************************************************************
CCCC ** CALCULATES THE CURRENT IN A DIODE OR P-N JUNCTION **
CCCC ** FOR A SERIES OF VALUES OF Vo **
CCCC ** SOLUTION IS FOUND BY SOLVING THE TRANSCENDENTAL EQUATION **
CCCC ** Vo=iR+(kT/e)ln(i/i(s)+1) **
CCCC ** where i(s)=reverse bias saturation current **
CCCC ** Vo =supply voltage **
CCCC ******************************************************************
DOUBLE PRECISION F
EXTERNAL F, C05ADF
CCCC Initialise Variables
DOUBLE PRECISION C,V,R,K,T,E,IS,TOL,X1,X2,X3,X4,ETA
INTEGER MNIT,J,IFAIL
C BLOCK DATA
C COMMON /PARAM/V,R,C,IS
C DATA R,K,T,E,IS,MNIT,TOL,ETA/1000,1.38E-23,300,1.6E-19,1.8E-15
C 1 ,50,1D-4,0.0/
C END
COMMON /PARAM/V,R,C,IS
R=1000
K=1.38D-23
T=300
E=1.6D-19
IS=1.8D-15
MNIT=50
TOL=1D-4
ETA=0.0
IFAIL=0
CCCC Print table
PRINT100,'VOLTAGE','CURRENT NAG ','CURRENT OWN'
DO 10 V=1,10
C=K*T/E
X1=0
X2=V/R
CALL SOLV(MNIT,X1,X2,J,X4,TOL)
C=K*T/E
X1=0
X2=V/R
CALL C05ADF (X1,X2,TOL,ETA,F,X3,IFAIL)
IF (IFAIL.NE.0) STOP
PRINT101,V,X3,X4
10 CONTINUE
100 FORMAT (3X,A,8X,A,5X,A,/)
101 FORMAT (3X,F5.2,7X,D12.5,9X,D12.5)
END
CCCC Solution of transcendental equations
CCCC by bisection method
SUBROUTINE SOLV(MNIT,X1,X2,J,X3,TOL)
DOUBLE PRECISION V,R,C,IS,X1,X2,X3,OLD,TOL
INTEGER J,MNIT
COMMON /PARAM/V,R,C,IS
DO 200 J=1,MNIT
X3=0.5*(X1+X2)
IF (F(X1)*F(X3).LT.0.0) THEN
X2=X3
ELSE
X1=X3
ENDIF
IF (J.NE.0) THEN
PRINT *,'OLD-F3 ',ABS(F(OLD)-F(X3))
C PRINT *,TOL
C PRINT *,'F3 ',F(X3)
IF (ABS(F(OLD)-F(X3)).LE.TOL) RETURN
ENDIF
OLD=X3
200 CONTINUE
END
CCCC FUNCTION
DOUBLE PRECISION FUNCTION F(X)
DOUBLE PRECISION X,V,R,C,IS
COMMON /PARAM/V,R,C,IS
F=X*R+C*LOG(X/IS+1.0)-V
END
ahh, uppercase programming, seems so retro these days. I can vaguely remember seeing the oscillations of a drumskin on the screen and I guess the following is the program that did it – though it used prospero’s own graphics routines so I can’t really check
CCCC ******************************************************************
CCCC ** DISPLAYS THE VIBRATIONS ON A DRUMSKIN AFTER BEING STRUCK **
CCCC ******************************************************************
INCLUDE 'F77PC'
REAL R,URT,A,K,T
INTEGER N
DIMENSION A(7),K(7)
PARAMETER (PI=3.14159)
DATA (A(I),I=1,6),(K(I),I=1,6)/0.5267,0.7074,0.7925,0.8433,
1 0.8642,0.8746,2.404,5.52,8.65,11.79,14.92,18.16/
CALL INISCR
DO 20 T=0,PI,0.1
DO 15 R=0,1.05,0.05
DO 10 N=1,7
URT=URT+A(N)*BESS(K(N)*R,1E-5)*COS(K(N)*T)
10 CONTINUE
C PRINT*,R,URT
CALL DISP(R,URT)
URT=0
15 CONTINUE
20 CONTINUE
END
REAL FUNCTION BESS(X,TOL)
REAL X,TOL,OLD,TERM,J
INTEGER N
TERM=1
J=1
N=0
IF (X.EQ.0.0) THEN
J=1
GOTO 120
ENDIF
DO WHILE (ABS(TERM-OLD).GE.TOL)
OLD=TERM
N=N+1
TERM=-(OLD*X*X)/(4.0*N*N)
J=J+TERM
END DO
120 BESS=J
END
SUBROUTINE INISCR
CALL INITSCREEN
CALL SETSCREENMODE(16)
CALL CLRVIDEO
END
February 14, 1991
BSc Fortran Numerics Course – Solving Equations
We did some very basic differential equation solving, simultaneous equations and curve fitting…
CCCC ******************************************************************
CCCC ** SOLUTION OF FIRST ORDER DIFFERENTIAL EQUATION **
CCCC ** USING EULERS METHODS **
CCCC ******************************************************************
EXTERNAL FUNC
REAL XO,XOUT,Y,H,FUNC
DATA XO,XOUT,Y/0.0,3.0,1.0/
PRINT *,'ENTER STEP SIZE '
READ *,H
CALL EULER (XO,XOUT,H,Y,FUNC)
PRINT *,'Y VALUE AT X=3.0 IS ',Y
END
REAL FUNCTION FUNC (X,Y)
REAL X,Y
FUNC=-Y
END
SUBROUTINE EULER (XO,XOUT,H,Y,FUNC)
REAL I,XO,XOUT,H,Y
DO 100 I=XO,XOUT,H
Y=Y+H*FUNC (XO,Y)
100 CONTINUE
END
and
CCCC ******************************************************************
CCCC ** SOLUTION OF FIRST ORDER DIFFERENTIAL EQUATION **
CCCC ** USING RUNGE-KUTTA METHODS **
CCCC ******************************************************************
EXTERNAL FUNC
REAL XO,XOUT,Y,H,FUNC
DATA XO,XOUT,Y/0.0,3.0,1.0/
PRINT *,'ENTER STEP SIZE '
READ *,H
CALL RK4 (XO,XOUT,H,Y,FUNC)
PRINT *,'Y VALUE AT X=3.0 IS ',Y
END
REAL FUNCTION FUNC (X,Y)
REAL X,Y
FUNC=-Y
END
SUBROUTINE RK4 (XO,XOUT,H,Y,FUNC)
REAL X,XO,XOUT,H,Y,C1,C2,C3,C4
DO 100 X=XO,XOUT,H
C1=H*FUNC (X,Y)
C2=H*FUNC (X+0.5*H,Y+0.5*C1)
C3=H*FUNC (X+0.5*H,Y+0.5*C2)
C4=H*FUNC (X+H,Y+C3)
Y=Y+(C1+2*C2+2*C3+C4)/6
100 CONTINUE
END
Simultaneous Equations
CCCC ******************************************************************
CCCC ** SOLUTION OF SIMULTANEOUS EQUATIONS **
CCCC ** BY GAUSSIAN ELIMINATION METHOD **
CCCC ** WITH PARTIAL PIVOTAL CONDENSATION **
CCCC ** AND CHECK FOR SINGULARITY **
CCCC ******************************************************************
CCCC Initialise variables
INTEGER I,J,K,N,L
REAL A,B,X,M,F,TOL,C,D,S
DIMENSION A(20,20),X(20),B(20)
COMMON A,B,X
CCCC Enter equations
PRINT*,'ENTER NO. OF EQUATIONS AND TOLERANCE '
READ*,N,TOL
DO 10 I=1,N
PRINT114,'ENTER COEFFICIENTS OF EQU.',I
READ*,(A(I,J),J=1,N),B(I)
10 CONTINUE
PRINT113,'EQUATIONS ENTERED'
CALL PR(N)
CCCC PARTIAL PIVOTAL CONDENSATION
DO 5 J=1,N
L=J
DO 2 I=1,N
IF (A(I,J).GT.A(L,J)) L=I
2 CONTINUE
DO 3 K=1,N
C=A(J,K)
A(J,K)=A(L,K)
A(L,K)=C
3 CONTINUE
D=B(J)
B(J)=B(L)
B(L)=D
5 CONTINUE
PRINT113,'PARTIAL PIVOTAL CONDENSATION'
CALL PR(N)
CCCC Gaussian elimination method
DO 50 K=1,N-1
DO 40 I=K+1,N
IF (ABS(A(K,K)).LT.TOL) STOP
M=A(I,K)/A(K,K)
S=0
DO 30 J=K,N
A(I,J)=A(I,J)-M*A(K,J)
S=S+A(I,J)**2
30 CONTINUE
B(I)=B(I)-M*B(K)
IF (S.LT.TOL) THEN
PRINT113,'SINGULARITY ERROR'
STOP
ENDIF
40 CONTINUE
50 CONTINUE
PRINT113,'TRIANGULAR FORM BEFORE BACK SUBSTITUTION'
CALL PR(N)
CCCC Back substitution & output of results
PRINT113,'SOLUTIONS'
X(N)=B(N)/A(N,N)
PRINT115,N,X(N)
DO 110 J=N-1,1,-1
DO 100 L=J+1,N
F=F+A(J,L)*X(L)
100 CONTINUE
X(J)=1/A(J,J)*(B(J)-F)
PRINT115,J,X(J)
F=0
110 CONTINUE
113 FORMAT (/,1X,A)
114 FORMAT (1X,A,I1,6X)
115 FORMAT (1X,'X(',I1,') = ',F9.5)
END
CCCC Display equations
SUBROUTINE PR(N)
INTEGER I,J,N
COMMON A(20,20),B(20),X(20)
DO 210 I=1,N
DO 200 J=1,N
PRINT215,A(I,J)
200 CONTINUE
PRINT216,B(I)
210 CONTINUE
215 FORMAT (1X,F9.4,4X,\)
216 FORMAT (1X,'=',2X,F9.4)
END
We were obviously using the Salford FTN77 compiler since there is a log file
SALFORD UNIVERSITY FTN77/386 VER. 1.67 A:\GAUSSPC.FOR 15:29:42 Tuesday, 22 January 1991
COMPILER OPTIONS: LISTING ANSI INTL NODCLVAR NOMAP NOCHECK LOGS DYNM OFFSET LGO NOPAGETHROW NOSILENT
NO_OPTIMISE NOWEITEK
I can remember having to go to the machine room in the Newton building to get the listings from the old 132 column printer (everywhere else was just epson 850’s which couldn’t print that wide.)
Finite Difference
CCCC ******************************************************************
CCCC ** CALCULATES THE TEMPERATURE IN A BLOCK OF HEATED MATERIAL **
CCCC ** USES INITIAL GUESSES AND RELAXATION TECHNIQUE APPLIED TO **
CCCC ** THE FINITE DIFFERENCE EQUATIONS AND DRAWS A GRAPH OF THE **
CCCC ** NUMBER OF ITERATIONS AGAINST THE RELAXATION COEFF. ALPHA **
CCCC ******************************************************************
CCCC Block is heated ...
CCCC 500 500 500 500
CCCC 100 ??? ??? 100
CCCC 100 ??? ??? 100
CCCC 100 100 100 100
REAL T,RES,RESMX,DECT,W,WMAX,WMIN,R,TT
INTEGER I,J,M,N,TMMIN,TMMAX,TNMIN,TNMAX,IT,PTS,NUM
DIMENSION T(10,10),TT(10,10),RES(10,10),PTS(50),W(50)
DATA ((TT(I,J),I=1,4),J=1,4),TMMIN,TMMAX,TNMIN,TNMAX/
1 4*500.0,100.0,2*300.0,2*100.0,2*200.0,5*100.0,1,4,1,4/
PRINT *,'ENTER RELAXATION COEFFICIENT BOUNDARIES '
READ *,WMIN,WMAX
DO 250 R=WMIN,WMAX,(WMAX-WMIN)/40.0
DO 5 I=1,4
DO 4 J=1,4
T(I,J)=TT(I,J)
4 CONTINUE
5 CONTINUE
DO 200 IT=0,200
CCCC FIND RESIDUALS
RESMX=0.0
DO 20 M=TMMIN+1,TMMAX-1
DO 10 N=TNMIN+1,TNMAX-1
RES (M,N)=T(M+1,N)+T(M-1,N)+T(M,N+1)+T(M,N-1)-4.0*T(M,N)
IF (ABS(RES(M,N)).GT.RESMX) RESMX=RES(M,N)
10 CONTINUE
20 CONTINUE
CCCC ALTER TEMPERATURE OF ALL POINTS BY RESMX/4.0*R
DECT=RESMX/4.0*R
DO 40 M=TMMIN+1,TMMAX-1
DO 30 N=TNMIN+1,TNMAX-1
T(M,N)=T(M,N)+DECT
30 CONTINUE
40 CONTINUE
IF (ABS(RESMX).LT.0.1 .AND. IT.GT.0) GOTO 230
200 CONTINUE
C 230 PRINT *,'REL COEFF',R,' STEPS ',IT
230 NUM=NUM+1
W(NUM)=R
PTS(NUM)=IT
250 CONTINUE
300 FORMAT (1X,I3,1X,12(F6.1))
CALL SCRSET
CALL PPOINT (NUM,W,PTS)
END
CCCC Initialise screen and clear it.
SUBROUTINE SCRSET
INCLUDE 'F77PC'
CALL INITSCREEN
CALL SETSCREENMODE (EGAGRAPH)
CALL CLRVIDEO
END
SUBROUTINE PPOINT (NPTS,X,Y)
INCLUDE 'F77PC'
INTEGER NPTS,Y,N,I,XP,YP
REAL X,XMAX,YMAX,SCALEX,SCALEY,XMIN,YMIN
DIMENSION Y(50),X(50)
CCCC Draw axis
CALL DRAW (20,330,620,330,WHITE)
CALL DRAW (20,330,20,10,WHITE)
CCCC Mark points
XMIN=1000
YMIN=1000
DO 390 N=1,NPTS
IF (X(N).GT.XMAX) XMAX=X(N)
IF (X(N).LT.XMIN) XMIN=X(N)
IF (Y(N).GT.YMAX) YMAX=Y(N)
IF (Y(N).LT.YMIN) YMIN=Y(N)
390 CONTINUE
SCALEX=600.0/(XMAX-XMIN)
SCALEY=320.0/(YMAX-YMIN)
DO 400,I=1,NPTS
XP=INT((X(I)-XMIN)*SCALEX+20)
YP=INT(330-(Y(I)-YMIN)*SCALEY)-2
CALL DRAW (XP-2,YP,XP+2,YP,WHITE)
CALL DRAW (XP,YP-2,XP,YP+2,WHITE)
400 CONTINUE
END
Least Squares Fitting
CCCC ******************************************************************
CCCC ** DRAWS A GRAPH OF DATA HELD IN FILE 'GDATA.DAT' **
CCCC ** AND DRAW LEAST SQUARES FIT **
CCCC ******************************************************************
CCCC Data file includes the number of points
CCCC and X,Y coordinates
INCLUDE 'F77PC'
INTEGER NPTS,I
REAL X,Y,GRAD,INTCPT
DIMENSION X(30),Y(30)
OPEN (5,FILE='A:\LSQ\GDATA.DAT',STATUS='OLD')
READ (5,*) NPTS
DO 10 I=1,NPTS
READ (5,*) X(I),Y(I)
10 CONTINUE
CALL LSFIT (X,Y,NPTS,GRAD,INTCPT)
CALL SCRSET
CALL GRAFIT (X,Y,NPTS,GRAD,INTCPT)
PRINT *,'GRADIENT=',GRAD
PRINT *,'INTERCEPT=',INTCPT
CLOSE (5)
END
CCCC Calculation of least squares fit to data
CCCC with coordinates held in array X,Y.
CCCC Returns GRAD-Gradient and INTCPT the intercept
CCCC on the y -axis
SUBROUTINE LSFIT (X,Y,NPTS,GRAD,INTCPT)
INTEGER NPTS,I
REAL X,Y,SX,SXX,SXY,SY,DEL,GRAD,INTCPT
DIMENSION X(*),Y(*)
DO 250 I=1,NPTS
SX=SX+X(I)
SXX=SXX+X(I)*X(I)
SXY=SXY+X(I)*Y(I)
SY=SY+Y(I)
250 CONTINUE
DEL=NPTS*SXX-(SX*SX)
GRAD=(NPTS*SXY-SX*SY)/DEL
INTCPT=(SXX*SY-SXY*SX)/DEL
END
CCCC Initialise screen and clear it.
SUBROUTINE SCRSET
INCLUDE 'F77PC'
CALL INITSCREEN
CALL SETSCREENMODE (EGAGRAPH)
CALL CLRVIDEO
END
CCCC Plots the data, draws the axis
CCCC and the line of best fit
SUBROUTINE GRAFIT (X,Y,NPTS,GRAD,INTCPT)
INCLUDE 'F77PC'
INTEGER I,NPTS,XP,YP,Y1,Y2
REAL X,Y,SCALEX,SCALEY,GRAD,INTCPT
DIMENSION X(*),Y(*)
CCCC Draw axis
CALL DRAW (20,330,620,330,WHITE)
CALL DRAW (20,330,20,10,WHITE)
CCCC Mark points
SCALEX=600/X(NPTS)
SCALEY=320/Y(NPTS)
DO 300,I=1,NPTS
XP=INT(X(I)*SCALEX+20)
YP=INT(330-Y(I)*SCALEY)
CALL DRAW (XP-2,YP,XP+2,YP,WHITE)
CALL DRAW (XP,YP-2,XP,YP+2,WHITE)
300 CONTINUE
CCCC Draw best fit line
Y1=INT(330-INTCPT*SCALEY)
Y2=INT(330-(GRAD*X(NPTS)+INTCPT)*SCALEY)
CALL DRAW (20,Y1,620,Y2,WHITE)
END
February 13, 1991
BSc Fortran Numerics Course – Fitting
I remember well this next program, I used it to fit a high order polynomial to some simple data – the result was a very wiggly line that ran through all the points! Though, I *think* this program is wrong – the original was in Propero and this looks like a half finished port of that program.
CCCC ******************************************************************
CCCC ** DRAWS A GRAPH OF DATA HELD IN FILE 'PDATA.DAT' **
CCCC ** AND DRAW LEAST SQUARES FIT TO POLYNOMIAL **
CCCC ******************************************************************
CCCC Data file includes the number of points
CCCC and X,Y coordinates
INTEGER DEG
DIMENSION XCOOR(20),YCOOR(20),COEFF(0:20,0:20),CONST(0:20)
OPEN (5,FILE='PDATA.DAT',STATUS='OLD')
READ (5,*) NPTS
DO 10 I=1,NPTS
READ (5,*) XCOOR(I),YCOOR(I)
10 CONTINUE
PRINT *,'ENTER DEGREE OF POLYNOMIAL '
READ *,DEG
CALL FITQ (XCOOR,YCOOR,NPTS,DEG,COEFF,CONST)
CALL GAUSPC(A,B,DEG+1,1E-3,G)
CALL SCRSET
CALL GRAFIT (X,Y,NPTS,DEG,MX,MY,MNX,MNY,A,B,G)
CALL GET_KEY@(K)
CALL TEXT_MODE@
CLOSE (5)
END
SUBROUTINE FITEQNS(XCOOR,YCOOR,NPTS,DEG,
* MATRIX,YPXSUM)
CCCC ==================================================================
CCCC Creates a set of simultaneous equations which when solved give the
CCCC best fit of a polynomial y = a + a(1) x + a(2) x**2 + ...
CCCC of degree DEG to the coordinates (XCOOR,YCOOR).
CCCC Returned is a matrix of the coefficients of the equations -
CCCC MATRIX(I,J) and the associated constants - YPXSUM(J).
CCCC Where I and J take values 0 -> DEG.
CCCC ==================================================================
CCCC DESCRIPTION OF VARIABLES
CCCC ------------------------
CCCC INTEGER DEG Degree of polynomial fit
CCCC NPTS Number of data points
CCCC REAL XCOOR() X - coordinates
CCCC YCOOR() Y - coordinates
CCCC MATRIX(,) Matrix of simultaneous eqn coeff
CCCC YPXSUM() Sum of Y times X to power n (n=0->DEG)
CCCC XPWSUM() Sum of X to power n (n=0->DEG)
CCCC ==================================================================
INTEGER DEG
REAL MATRIX
DIMENSION XCOOR(*),YCOOR(*),
* MATRIX(0:100,0:*),YPXSUM(0:*),XPWSUM(0:*)
C
C Find the sum of the x-coordinate to the power of n
C where n varies between 0 and 2 times the degree of polynomial
C for all coordinates.
C
XPWSUM(0)=NPTS
DO 20 I=1,2*DEG
DO 10 J=1,NPTS
XPWSUM(I)=XPWSUM(I)+XCOOR(J)**I
10 CONTINUE
20 CONTINUE
C
C Find the sum of (x-coordinate to the power n) times by the
C y-coordinate for all coordinates.
C
DO 110 I=0,DEG
DO 100 J=1,NPTS
YPXSUM(I)=YPXSUM(I)+YCOOR(J)*XCOOR(J)**I
100 CONTINUE
110 CONTINUE
C
C Create matrix of coefficients of simultaneous equations
C of the form a(0,0) + a(1,0) + ... + a(i,0) = const(0)
C a(0,1) + a(1,1) + ... + a(i,1) = const(1)
C : : : :
C a(0,j) + a(1,j) + ... + a(i,j) = const(j)
C
C Note the const() is YPXSUM() and is already evaluated.
C
DO 210 J=0,DEG
DO 200 I=0,DEG
MATRIX(I,J)=XPWSUM(I+J)
200 CONTINUE
210 CONTINUE
END
CCCC Solve the best fit equations
SUBROUTINE GAUSPC (A,B,N,TOL,G)
INTEGER I,J,K,N,L
REAL A,B,X,M,F,TOL,C,D,G
DIMENSION A(20,20),G(20),B(20)
CCCC PARTIAL PIVOTAL CONDENSATION
DO 320 J=1,N
L=J
DO 300 I=1,N
IF (A(I,J).GT.A(L,J)) L=I
300 CONTINUE
DO 310 K=1,N
C=A(J,K)
A(J,K)=A(L,K)
A(L,K)=C
310 CONTINUE
D=B(J)
B(J)=B(L)
B(L)=D
320 CONTINUE
CCCC Gaussian elimination method
DO 450 K=1,N-1
DO 400 I=K+1,N
IF (ABS(A(K,K)).LT.TOL) STOP
M=A(I,K)/A(K,K)
DO 350 J=K,N
A(I,J)=A(I,J)-M*A(K,J)
350 CONTINUE
B(I)=B(I)-M*B(K)
400 CONTINUE
450 CONTINUE
CCCC Back substitution
G(N)=B(N)/A(N,N)
DO 110 J=N-1,1,-1
f=0
DO 100 L=J+1,N
F=F+A(J,L)*G(L)
100 CONTINUE
G(J)=1/A(J,J)*(B(J)-F)
110 CONTINUE
END
CCCC Initialise screen and clear it.
SUBROUTINE SCRSET
CALL VGA@
CALL CLEAR_SCREEN@
END
CCCC Plots the data, draws the axis
CCCC and the line of best fit
SUBROUTINE GRAFIT (X,Y,NPTS,DEG,MX,MY,MNX,MNY,A,B,G)
INTEGER I,NPTS,XP,YP,Y1,Y2,OLDX,OLDY,DEG
REAL X,Y,SCALEX,SCALEY,MX,MY,MNX,MNY,A,B,G,C,D
DIMENSION X(20),Y(20),A(20,20),B(20),G(40)
IF (MNX.GT.0.0) MNX=0.0
IF (MNY.GT.0.0) MNY=0.0
SCALEX=600.0/(MX-MNX)
SCALEY=310.0/(MY-MNY)
CCCC Draw axis
XP=INT(20-MNX*SCALEX)
YP=INT(330+MNY*SCALEY)
CALL DRAW_LINE@(XP,20,XP,330,15)
CALL DRAW_LINE@(20,YP,620,YP,15)
CCCC Mark points
DO 500,I=1,NPTS
XP=INT(X(I)*SCALEX+20-MNX*SCALEX)
YP=INT(330+MNY*SCALEY-Y(I)*SCALEY)
CALL DRAW_LINE@(XP-2,YP,XP+2,YP,15)
CALL DRAW_LINE@(XP,YP-2,XP,YP+2,15)
500 CONTINUE
CCCC Draw Curve
DO 620 C=MNX,MX+0.005,0.01
D=0
DO 610 I=1,DEG+1
IF (C.EQ.0.0) THEN
D=G(1)
GOTO 610
ENDIF
D=D+G(I)*C**(I-1)
610 CONTINUE
XP=INT(C*SCALEX+20-MNX*SCALEX)
YP=INT(330+MNY*SCALEY-D*SCALEY)
IF (OLDX.EQ.0.0) GOTO 615
CALL DRAW_LINE@(OLDX,OLDY,XP,YP,15)
615 OLDX=XP
OLDY=YP
620 CONTINUE
PRINT *,'DEGREE OF POLYNOMIAL FIT ',DEG
END
Iterative Root Finding
CCCC ******************************************************************
CCCC ** EVALUATION OF SQUARE ROOT **
CCCC ** BY ITERATIVE METHOD **
CCCC ** **
CCCC ** Improved guess= 0.5 (guess+number/guess) **
CCCC ** **
CCCC ******************************************************************
CCCC Initialise variables
REAL NUM,GU,TOL,VAL
INTEGER FLAG,MNIT,N
CCCC Enter data
PRINT*,'ENTER NUMBER, GUESS, TOLERANCE, MAX.ITERATIONS '
READ*,NUM,GU,TOL,MNIT
CCCC Call square root routine
CALL SQ(NUM,TOL,MNIT,GU,FLAG,VAL,N)
IF (FLAG.GE.2) THEN
IF (FLAG.EQ.2) PRINT *,'NUMBER TOO SMALL OR NEGATIVE'
IF (FLAG.EQ.3) PRINT *,'GUESS TOO SMALL OR NEGATIVE'
STOP
ELSE
IF (FLAG.EQ.1) PRINT *,'WARNING:MAXIMUM NO. OF ITERATIONS',
1 ' WAS REACHED'
ENDIF
PRINT*,'VALUE IS ',VAL,' NO. OF ITERATIONS',N
END
CCCC Square root routine
SUBROUTINE SQ(NUM,TOL,MNIT,GU,FLAG,VAL,N)
REAL NUM,TOL,GU,VAL
INTEGER MNIT,FLAG,N
IF (NUM.LT.TOL) THEN
FLAG=2
RETURN
ELSE
IF (GU.LT.TOL) THEN
FLAG=3
RETURN
ENDIF
ENDIF
DO 10 N=0,MNIT
VAL=(GU+NUM/GU)*0.5
IF (ABS(GU-VAL).LT.TOL) RETURN
GU=VAL
10 CONTINUE
FLAG=1
RETURN
END
But I think writing our own routines was just an exercise; we used the NAG library too
cccc*******************************************************************
cccc Calculates the best fit of data in requested file
cccc to equation given in lsfun1, uses NAG E04FDF
cccc*******************************************************************
cccc===================================================================
c Equation is the number of counts from parent-daughter decay
c Find best fit of data using constants lambda a & b the decay
c rates and c0 the initial decay rate.
c filenm - filename holding experimental data
c t() - time
c ct() - counts at time t
c m - number of data points
c n - number of unknown constants in equation
cccc===================================================================
integer t(50), ct(50)
parameter (m=41, n=3, liw=1, lw=7*n+n*n+2*m*n+3*m+n*(n-1)/2)
double precision cttmp(50), iw(liw), w(lw), x(n)
character filenm*8
external e04fdf
common t, ct
cccc get filename for experimental data
print *,'Enter filename containing experimental data '
read (*,1000) filenm
open (5,file=filenm)
cccc read in data file and convert counts to integer
do 10 i=1,m
read (5,*) t(i), cttmp(i)
ct(i)=int(cttmp(i))
10 continue
cccc initial guesses
x(1)=1.0d-3
x(2)=3.0d-3
x(3)=10040.0d0
ifail=1
call e04fdf(m, n, x, fsumsq, iw, liw, w, lw, ifail)
print *,'Minimum point is ',(x(j),j=1,n)
c print 1001,'Minimum point is ',x(1),x(2),x(3)
close (5)
1000 format (a)
1001 format (1x,'lambda a ',d9.3,' lambda b ',d9.3,
* ' initial count ', d9.3)
end
subroutine lsfun1(m, n, xc, fvecc)
cccc*******************************************************************
cccc Calculates the difference between the predicted count rate
cccc using the estimates constants lambda a,b and the initial count
cccc*******************************************************************
cccc===================================================================
c la - lambda a
c lb - lambda b
c c0 - initial count
c fvecc - difference between calculated and experimental values
cccc===================================================================
integer t, ct
double precision fvecc(50), xc(3), la, lb, c0
dimension t(50),ct(50)
common t, ct
la=xc(1)
lb=xc(2)
c0=xc(3)
do 10 i=1,m
fvecc(i)=((lb*2-la)*exp(-la*t(i))-
* lb*exp(-lb*t(i)))*c0/(lb-la)-ct(i)
10 continue
end
Interpolation, Extrapolation and Bisection
CCCC ******************************************************************
CCCC ** PROGRAM TO DEMONSTRATE A VARIETY OF METHODS FOR SOLVING **
CCCC ** TRANSCENDENTAL EQUATIONS **
CCCC ******************************************************************
CCCC
CCCC THIS PROGRAM GIVES THE ANGLE OF A GALVONOMETER (THI) FOR A GIVEN
CCCC CURRENT THROUGH IT BY SOLVING THE TRANSCENDENTAL EQUATION
CCCC THI=(NAB/C)ICOS(THI)
CCCC THIS PROGRAM HOLDS DATA FOR THE VARIABLES BUT MAY BE USED WITH
CCCC OTHER VALUES BY REMOVING THE COMMENT MARKS IN THE MARGIN BELOW
CCCC AND THE DATA STATEMENT.
CCCC Initialise variables
REAL C,A,B,TOL,D,THI,X1,X2,I,PIB2,E,THI1,THI2,THI3,THI4
INTEGER N,STEP,MNIT,J,J1,J2,J3,J4
PARAMETER (PIB2=3.141592/2)
CCCC Enter variable values
C PRINT *,'ENTER C,N,A,B,LENGTH OF TABLE,MAX.NO ITERATIONS AND '
C 1 ,'TOLERANCE'
C READ*,C,N,A,B,STEP,MNIT,TOL
DATA C,N,A,B,STEP,MNIT,TOL/5,200,0.15,0.05,5,200,0.001/
D=C/(N*A*B)
CCCC Print table heading
PRINT99
PRINT100
PRINT101
PRINT 99
CCCC Submit current values to subroutines
CCCC and display solutions
DO 10 I=0,4*D,4*D/STEP
E=I/D
X1=0
X2=PIB2
IF (F(X1,E)*F(X2,E).GT.0.0) THEN
PRINT *,I,'NO ROOTS'
GOTO 10
ENDIF
CALL BISECT (E,MNIT,X1,X2,J,THI)
J1=J
THI1=THI
X1=0
X2=PIB2
CALL NEWT (E,MNIT,X1,J,THI,TOL)
J2=J
THI2=THI
X1=0
X2=PIB2
CALL INTERP (E,MNIT,X1,X2,J,THI,TOL)
J3=J
THI3=THI
X1=0
X2=PIB2
CALL EXTRAP(E,MNIT,X1,X2,J,THI,TOL)
J4=J
THI4=THI
PRINT102,I,J1,THI1,J2,THI2,J3,THI3,J4,THI4
10 CONTINUE
PRINT99
99 FORMAT (1X,62('-'))
100 FORMAT (1X,'CURRENT',3X,'BISECTION',5X,'NEWTON',3X,
1 'INTERPOLATION',3X,'EXTRAPOLATION')
101 FORMAT (11X,'N',3X,'THI',6X,'N',2X,'THI',7X,'N',4X,'THI',
1 8X,'N',3X,'THI')
102 FORMAT (1X,F7.4,1X,I3,1X,F7.4,2X,I3,1X,F7.4,2X,I3,1X,F7.4,
1 5X,I3,1X,F7.4)
END
CCCC Solution by bisection method
SUBROUTINE BISECT(E,MNIT,X1,X2,J,X3)
REAL D,X1,X2,X3,E,OLD,TOL
INTEGER MNIT,J
DO 200 J=0,MNIT
X3=0.5*(X1+X2)
IF (F(X1,E)*F(X3,E).LT.0.0) THEN
X2=X3
ELSE
X1=X3
ENDIF
IF (J.NE.0) THEN
IF (ABS(F(OLD,E)-F(X3,E)).LE.TOL) RETURN
ENDIF
OLD=X3
200 CONTINUE
END
CCCC Solution by Newton's method
SUBROUTINE NEWT(E,MNIT,X1,J,X2,TOL)
REAL E,X1,X2,GRAD,TOL
INTEGER J,MNIT
DO 300 J=0,MNIT
GRAD=(F(X1+0.01,E)-F(X1,E))/0.01
X2=X1-(F(X1,E)/GRAD)
IF (J.NE.0) THEN
IF (ABS(F(X2,E)-F(X1,E)).LE.TOL) RETURN
ENDIF
X1=X2
300 CONTINUE
END
CCCC Solution by interpolation
SUBROUTINE INTERP(E,MNIT,X1,X2,J,X3,TOL)
REAL X1,X2,X3,E,GRAD,OLD,TOL
INTEGER J,MNIT
DO 400 J=0,MNIT
GRAD=(F(X2,E)-F(X1,E))/(X2-X1)
X3=X1-F(X1,E)/GRAD
IF (F(X1,E)*F(X3,E).LT.0.0) THEN
X2=X3
ELSE
X1=X3
ENDIF
IF (J.NE.0) THEN
IF (ABS(F(OLD,E)-F(X3,E)).LE.TOL) RETURN
ENDIF
OLD=X3
400 CONTINUE
END
CCCC Solution by extrapolation
SUBROUTINE EXTRAP(E,MNIT,X1,X2,J,X3,TOL)
REAL X1,X2,X3,E,GRAD,OLD,TOL
INTEGER J,MNIT
DO 500 J=0,MNIT
IF ((X2-X1).EQ.0.0) THEN
J=0
X3=0
RETURN
ENDIF
GRAD=(F(X2,E)-F(X1,E))/(X2-X1)
X3=X1-F(X1,E)/GRAD
IF (ABS(X3-X1).LT.ABS(X3-X2)) THEN
X2=X3
ELSE
X1=X3
ENDIF
IF (J.NE.0) THEN
IF (ABS(F(OLD,E)-F(X3,E)).LE.TOL) RETURN
ENDIF
OLD=X3
500 CONTINUE
END
CCCC Form of transcendental equation
FUNCTION F(X,E)
REAL F,X,E
F=X-E*COS(X)
END
There were perhaps another dozen or so other programs for e.g. calculating a certain branch of arctan by a formula etc. but not very interesting by today’s coding standards.
