* Move particles (fish) according to external force field (current) * * Inputs : * nfish = number of fish * fishp = array of initial fish positions (stored as complex numbers) * fishv = array of initial fish velocities (stored as complex numbers) * fishm = array masses of fish * tfinal = final time for simulation * zoom = range of display field * * Outputs: * time intervals and root mean square velocity at each step * plot of fish positions after each time step * * Algorithm: integrate using Euler's method * program fish1 include '/usr/cm5/cmf/include/cm/timer-fort.h' include '/usr/cm5/cmf/include/cm/CMF_defs.h' * Define the size of the display window integer, parameter :: m = 256 integer, parameter :: disp_interval = 10 include 'color_display.inc' * Max. No. of fish integer, parameter :: nfish = 10000 complex fishp(nfish), fishv(nfish), force(nfish), accel(nfish) real fishm(nfish), fishspeed(nfish) CMF\$ layout fishp(:news),fishv(:news),fishm(:news) CMF\$ layout force(:news),accel(:news) integer step,i real dt, t, tfinal, mnsqvel, zoom complex, parameter :: ii = (0,1) C C Initialize fishes, time, and display parameters C force = (0.,0.) accel = (0.,0.) forall (i=1:nfish) fishp(i) = (float(i)*2.0/nfish)-1.0 fishv = ii*fishp forall (i=1:nfish) fishm(i) = 1.0+float(i)/nfish tfinal = 30. zoom = 40. * * Loop over time steps step = 0 dt = .01 t = 0.0 do while (t < tfinal) t = t + dt * * Update fish positions and velocities fishp = fishp + dt*fishv call compute_current(force,fishp) accel = force/fishm fishv = fishv + dt*accel * Update time step fishspeed = abs(fishv) mnsqvel = sqrt(sum(fishspeed*fishspeed)/nfish) print 100,t,mnsqvel dt = .1*maxval(fishspeed) / maxval(abs(accel)) * * Show fish positions if (mod(step,disp_interval) .eq. 0) then call show_fish(fishp,nfish,zoom,m) endif step = step + 1 enddo 100 format(2f20.10) stop end subroutine compute_current(force,fishp) complex force(:),fishp(:) complex, parameter :: ii = (0,1) CMF\$ layout force(:news),fishp(:news) force = (3,0)*(fishp*ii)/(max(abs(fishp),0.01)) - fishp end subroutine show_fish(fishp,nfish,zoom,m) complex fishp(nfish) real zoom integer x(nfish), y(nfish), show(m, m) integer j * Scale the coordinates according to the zoom, and shift the origin x = INT(( real(fishp)+zoom)/(zoom)*(m/2)) y = INT((aimag(fishp)+zoom)/(zoom)*(m/2)) show = 0 forall(j=1:nfish) show(x(j),y(j)) = 1 * * enlarge each pixel for better visibility show(:,1:m-1) = show(:,1:m-1) + show(:,2:m) show(1:m-1,:) = show(1:m-1,:) + show(2:m,:) call display(show,m,m) end