Rcjp's Weblog

September 12, 2007


Filed under: physics — rcjp @ 4:09 pm

Everytime I look at this image, I have to stare at it for several minutes in sheer astonishment. Kudos to NASA for making these images available.

April 12, 2007

Circular Halo

Filed under: physics — rcjp @ 12:53 pm


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))

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))

much better, phew.

April 5, 2007


Filed under: physics, python — rcjp @ 3:56 pm

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
import gtk, math

class Pattern(object):
    x, y = 0, 0
    alpha = 76.11
    def __init__(self, xsize, ysize):
        self.xsize, self.ysize = xsize, ysize
        window = gtk.Window(gtk.WINDOW_TOPLEVEL)
        window.connect("destroy", lambda w: gtk.main_quit())
        vbox = gtk.VBox(homogeneous=False, spacing=5)
        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)

    def slider_changed(self, slider):
        self.alpha = slider.get_value()

    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)
            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:

    def area_expose_cb(self, area, event):
        blue = self.area.get_colormap().alloc_color("#0000FF")
        pointgc = self.area.window.new_gc()
        area.window.draw_point(pointgc, 20, 10)
        self.area.window.draw_points(pointgc, self.points)

def main():
    return 0

if __name__ == '__main__':
    Pattern(600, 500)

November 16, 2006

Pattern Formation

Filed under: physics — rcjp @ 10:09 am

cccc program from 'Models of Biological Pattern Formation' Hans Meinhardt

! 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

      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)
         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)
         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

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)
         print *


November 7, 2006

Series Solution Method in maxima

Filed under: physics — rcjp @ 5:56 pm

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

  and using series method
eq1: diff(q(t), t)= p(t)/m;
eq2: diff(p(t), t)= -m*w^2*q(t);
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 */


/* Note rearrange eqn=0 */

eqns:append(cons(subst(0,t,p) = 0,coeffs1),cons(subst(0,t,q) = 1,coeffs2))$

or inside an emaxima session:

% -*-EMaxima-*-
\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

  eq1: diff(q(t), t)= p(t)/m;
  eq2: diff(p(t), t)= -m*w^2*q(t);
  atvalue(p(t), t=0, 0)$
  atvalue(q(t), t=0, 1)$
  ans: desolve([eq1, eq2], [q(t), p(t)]);

To do that with series expansions only:

  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));


November 4, 2006


Filed under: physics — rcjp @ 4:21 pm

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-*-
    \usepackage{a4wide}  % bit more room!
    \usepackage{amsmath} % used to get e.g. \boxed{...}
    \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$$

      'integrate(cosh(a*x)*f(t-x), x, 0, t) + b*f(t)=t**3;
      res: linsolve([%],['laplace(f(t),t,s)]);

    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

      dpart(res[1], 2);
      ilt(part(res[1], 2),s,t);


and a slightly more complex example…

    % -*-EMaxima-*-
    %\usepackage{a4wide}  % bit more room!
    \usepackage{amsmath} % used to get e.g. \boxed{...}

    \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$

      eqn_init: 'diff(u,t,2)+u=epsilon*u^2;
      v: u0 + epsilon*u1 + epsilon^2*u2;
      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);

February 23, 2005


Filed under: physics — rcjp @ 12:04 pm


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

Filed under: physics — rcjp @ 8:31 pm


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

Filed under: physics — rcjp @ 11:57 am

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

      9  *   4  =  36
      7  *  10  =  70
     50  -  36  =  14
     70  *  14  = 980
    980  -   9  = 971

      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
    end module solutions

    program countdown
      use solutions

      numbers=(/50, 7, 5, 1, 8, 3/)

      numbers=(/50, 9, 9, 4, 7, 10/)
      soluop=' '
      call pair(numbers, ncards, solunum, solunum1, soluop, ncards, ntarget)
      print *, 'at least ',nsols-1, ' solutions'
    end program countdown

    integer function combine(n, m, op)
      integer::n, m

      select case (op)
      case ('*')
      case ('+')
      case ('-')
      case ('/')
         if (mod (n,m) /= 0) combine=0
      case default
      end select
      if (combinesolsize*2) then
            print *, '***SOLUTION***'
            do m=solsize, 2, -1
               select case (soluop(m))
               case ('*')
               case ('+')
               case ('-')
               case ('/')
               end select
               print '(i4, 2x, a, i4, 2x, a, i4)', n1,soluop(m), n2, '=', n3
    end subroutine printsolution

October 27, 1995

Teaching MSc Numerical Methods

Filed under: physics — rcjp @ 11:51 am

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…

Older Posts »

Blog at WordPress.com.