Rcjp's Weblog

November 16, 2006

Pattern Formation

Filed under: physics — rcjp @ 10:09 am

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

Pattern

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: