program posphere
implicit none
c----------------------------------------------------------------------
c This program generates the time progression of a system of n equal
c charges constrained to a sphere of radius one.
c What must be done in this program:
c 1) generate n x,y,z position vectors - see gennpart subroutine
c (gennpart.f)
c 2) adjust the moduli of the position vectors so they lie on the
c sphere. Note: This subroutine reflonsphere is used
c thereafter to ensure the particles are constrained to the
c sphere. (reflonsphere.f)
c 3) generate velocities (in the case of this program, all of
c the components of every particles initial velocity is set
c equal to one).
c 4) generate initial forces (from the initial positions).
c see force.f for details of the force calculations
c 5) calculate the new velocities, forces, and positions
c starting with the inital velocities, forces, and
c positions. To find the new values we use a FDA of
c the equation of motion F = ma. a is simply the second
c derivative of position and the forces, acting on the system,
c are the coulomb forces (between each particle) and a damping
c force (-gamma*v) necessary to bring the system to equilibrium.
c To solve this FDA we will split it up into two differential
c equations. The leap-frog method will be used. In the leap-
c frog method, we calculate velocities at half time steps and
c positions at full time steps. The two FD equations that were
c derived are:
c v(n+.5) = dt*(sum of forces) + ( 1 - dt*gamma/2)*v(n-.5) *
c
c & (1 + dt*gamma/2)^-1 (1)
c
c r(n+1) = r(n-1) + dt * v(n+.5) (2)
c
c All velocities and forces will be made tangent to the sphere
c to prevent the charges from flying away to infintiy. (see
c tangent.f subroutine for details)
c 6) create output. Several possibilities. The (get name of sub)
c subroutine, written by Dr. Matthew W. Choptuik animates the
c position vectors that his program can output. Otherwise we
c can print the epsilon values at given time steps. Epsilon
c is the magnitude of the difference vector of position vectors
c at subsequent times. A threshold epsilon is in the do while
c loop of this program. As the particles approach equilibrium,
c the difference vector approaches zero (as does epsilon) by
c changing the value of epsilon that terminates the program, we
c change the accuracy of the equilibrium values. The most
c interesting output are the final equilibrium partical distances
c with these we can build the configurations.
c
c Declared variables:
c
c Real*8 -->
c
c r8arg - External function (p329f) that fetches real*8 values
c from the command line. Its first argument is the place
c of the number to be fetched (its # on the command line). Its
c second argument is a default value.
c
c dt - The time mesh spacings. User defined as maxtime /
c nsteps.
c
c ddotprod - external function that returns the dot product of
c two vectors sent in the function call. (see ddotprod.f for
c details.
c
c radius - holds the radius of the inidividual charge spheres.
c For use in sphplot.c. This subroutine was instrumental
c in examining equilibrium configurations. It animates the
c time progression of the charges on the surface of the sphere.
c
c posn - 2d array that stores n(number of particles) length 3
c position vectors.
c
c vel - 2d array that stores the velocities corresponding
c to the above positions (3 x n).
c
c gamma - The damping term.
c
c eps - value of epsilon given by the user on the command
c line. Default value = 1.0d-16.
c
c netforce - stores the netforce on the particles from the
c other n-1 identically charged particles. (3 x n)
c
c diff - 3 x 1 work array. temperarily stores the differnce
c vectors (position), for use in calculation of the distances
c between particles at equilibrium.
c
c epsilon - explained above.
c
c step100 - used to output, to the user, what time step the
c program is on. (multiples of 100).
c
c ngd2p1 - constant for use in finding new velocities. Equals
c 1 - dt*gamma/2.
c
c gd2p1 - constant for use in finding new velocities. Equals
c 1 + dt*gamma/2.
c
c chk1 and chk2 - arrays used for storing position vectors of
c subsequent times. Used to calculate epsilon.
c
c Integers -->
c
c i4arg - integer equivalent of r8arg (above).
c
c iargc - returns the number of arguments on the command line.
c (p329f)
c
c i, j, k, and l - counter variables.
c
c maxn - maximum value of all 3 x n arrays is 3 x maxn.
c Adjustable array value.
c
c n - the number of particles.
c
c maxtime - the last value in the time mesh.
c nsteps - the number of time values in the mesh. (see dt
c above).
c
c Special thanks to Dr. Matthew W. Choptuik, for all of his
c advice and for the use of the dumpvec subroutine that made
c debugging easier. sphplot was also authored by Dr. Choptuik
c and it made examination of the equilibrium configurations
c (visually) a breeze.
c
c----------------------------------------------------------------------
integer i4arg, iargc
real*8 r8arg, dt, ddotprod, radius
integer i, j, k, l
integer maxn, n, maxtime, nsteps
parameter (maxn = 600 000)
real*8 posn(3,maxn), vel(3,maxn), gamma, eps
real*8 netforce(3,maxn), diff(3), epsilon, step100
real*8 ngd2p1, gd2p1, chk1(3), chk2(3)
c----------------------------------------------------------------------
c Get level, n, and gamma from the command line. Note level and gamma
c have default values.
c----------------------------------------------------------------------
if (iargc() .lt. 1 .or. iargc() .gt. 5) goto 900
n = i4arg(1,-1)
maxtime = i4arg(2,8)
nsteps = i4arg(3,20)
gamma = r8arg(4,2.0d0)
eps = r8arg(5,1.0d-16)
radius = 0.05d0
c----------------------------------------------------------------------
c Call gennpart to generate n position vectors. Then call reflonsphere
c to give them modulus 1.
c----------------------------------------------------------------------
write(0,*)eps
call gennpart(posn,n)
call reflonsphere(posn,n,1)
c----------------------------------------------------------------------
c Initialize dt (maxtime/nsteps) and the constants ngd2p1 and gd2p1.
c These are constants in the FDA of the new positions.
c----------------------------------------------------------------------
dt = maxtime * nsteps**(-1.0d0)
ngd2p1 = 1 - gamma * 0.5d0 * dt
gd2p1 = 1.0d0/(1 + gamma * 0.5d0 * dt)
c----------------------------------------------------------------------
c Generate random velocities and make them tangent to the surface of
c the sphere.
c----------------------------------------------------------------------
do j = 1, n
do i = 1,3
vel(i,j) = 1.0d0
end do
end do
do i = 1, n
call tangent(vel(1,i),posn(1,i),3)
end do
c----------------------------------------------------------------------
c Find initial force and make the force tangent to the sphere.
c----------------------------------------------------------------------
call force(posn,netforce,n)
do i = 1, n
call tangent(netforce(1,i),posn(1,i),3)
end do
c----------------------------------------------------------------------
c Now we are prepared to iterate the system nsteps times with the
c time of each step being dt. chk1 and chk2 are found below. they
c hold the position vector of particle one of subsequent times. i.e.
c chk2 holds r1(dt*n) and chk1 holds r(dt*(n-1)) where dt*n is the
c newly calculated position vector. There are two ways to use this
c program. One in which the user inputs the total # of steps to be
c carried out and one in which the program continually subtracts the
c vectors chk1 and chk2 and compares the magnitude of their difference.
c As magnitude(chk1-chk2) --> 0; epsilon --) 0. In the do while
c statement below we can control the accuracy of the equilibrium
c positions of the particles by changing epsilon. eps stores the user
c defined epsilon (default = 1.0d-16). The key parameter in this
c program is (dt*gamma)/2. As with all key paramters of FDA's the
c user must be careful about the values he/she inputs for this
c parameter in order for convergence to occur. (dt = maxtime/nsteps)
c (dt*gamma)/2 is the courant number of the approximation. It's value
c and the value of epsilon will determine the accuracy of the final
c configuration.
c----------------------------------------------------------------------
chk2(1) = 1
chk2(2) = 1
chk2(3) = 1
c do k = 1, nsteps
k = 2
epsilon = 1
do while (epsilon .gt. eps)
do i = 1,3
chk1(i) = chk2(i)
end do
c---------------------------------------------------------------------
c Calculate the new velocities using the finite difference approx..
c Note: dmmsmsa is again used because it need only the netforce
c and vel arrays and returns the new vel array. The new vel array
c is, of course, A linear combination of the sums of vel and netforce
c---------------------------------------------------------------------
call dmmsmsa(vel,netforce,vel,gd2p1*ngd2p1,dt*gd2p1,
& 3,n)
c---------------------------------------------------------------------
c Find new positions with the new velocities and reflect those new
c positions onto the sphere. Also, here I find the new value of
c chk2 and subtract chk2 from chk1. The difference is stored in
c chk1. epsilon (the magnitude of the difference) is also computed
c here.
c---------------------------------------------------------------------
call dmmsmsa(posn,vel,posn,1.0d0,dt,3,n)
call reflonsphere(posn,n,1)
call dmmsmsa(posn(1,1),posn(1,2),chk2,1.0d0,-1.0d0,3,1)
call dmmsmsa(chk1,chk2,chk1,1.0d0,-1.0d0,3,1)
epsilon = sqrt(ddotprod(chk1,chk1,3))
c---------------------------------------------------------------------
c Make the new velocities tangent to the new position vectors
c---------------------------------------------------------------------
do l = 1, n
call tangent(vel(1,l),posn(1,l),3)
end do
c---------------------------------------------------------------------
c Calculate new force (with new position vectors)
c---------------------------------------------------------------------
call force(posn,netforce,n)
c---------------------------------------------------------------------
c Make the new forces tangent to the sphere
c---------------------------------------------------------------------
do l = 1, n
call tangent(netforce(1,l),posn(1,l),3)
end do
c---------------------------------------------------------------------
c Here we output to the user, in increments of 100, the current time
c step. Here also is the end of the main loop.
c---------------------------------------------------------------------
step100 = dmod(k,100.0d0)
if (step100 .eq. 0.0d0) then
write(0,*)k,'th time step'
write(0,*)epsilon,' <<epsilon>>'
end if
k = k + 1
c---------------------------------------------------------------------
c Call the subroutine sphplot when we want to see the system come
c to equilibrium visually.
c---------------------------------------------------------------------
c call sphplot(posn,n,radius)
end do
c---------------------------------------------------------------------
c Call the subroutine sphplot after equilibrium has occured and with
c the number of charges sent with a negative value. This creates a
c tumbling effect and makes it easier to see the configurations
c 3 - dimensionally.
c---------------------------------------------------------------------
c call sphplot(posn,-n,radius)
c---------------------------------------------------------------------
c Here we output the final position, the magnitude of each particle's
c velocity and the distances between the particles in the final
c configuration.
c---------------------------------------------------------------------
write(0,*)
write(0,*)
write(0,*) '<<final positions:>>'
do i = 1, n
write(0,*)posn(1,i),posn(2,i),posn(3,i)
end do
write(0,*)
write(0,*)
do i = 1, n
write(0,*)i,'th particles velocity (magnitude)'//': ',
& sqrt(ddotprod(vel(1,i), vel(1,i),3))
end do
write(0,*)
write(0,*)
do i = 1, n
do j = i, n
if (i .ne. j) then
call dmmsmsa(posn(1,i),posn(1,j),
& diff,1.0d0,-1.0d0,3,1)
write(0,*)'**********************'
write(0,*)'distance between particles:',i,' and',j,':'
write(0,*) sqrt(ddotprod(diff,diff,3))
write(0,*)'**********************'
end if
end do
end do
c---------------------------------------------------------------------
c Normal (no error) exit
c---------------------------------------------------------------------
stop
c----------------------------------------------------------------------
c Program comes here if a usage error occured
c Print usage error and exit
c----------------------------------------------------------------------
900 continue
write(0,*) 'Usage: posphere <n> [<maxtime> <nsteps> <gamma> '//
& '<eps>]'
write(0,*) 'Default values: 8 20 2.0d0 '//
& '1.0d-16'
stop
end
c**********************************************************************
c**********************************************************************
c**********************************************************************
subroutine force(r,f,n)
implicit none
c---------------------------------------------------------------------
c force:
c This subroutine calculates the coulomb force between two
c particles of identical charge. Note: the force is calculated
c in electrostatic units and the mass and charge of the particles
c has been arbitrarily set equal to one. The force on each charge
c (sum of the individual charges) is calculated and placed into
c f. This force array is returned to the main program.
c
c Variables:
c
c real*8 -->
c
c r(3,n) --> 3 x n array containing the current position vectors
c of n charges. Sent from the main program.
c
c f(3,n) --> 3 x n array in which the sum of the forces on
c each particles is calculated. The column of the array
c corresponds to the total force on the nth particle.
c
c rij(3) --> column vector where the vector from particle j to
c particle i is calculated. It is used in finding individual
c forces between particles i and j.
c
c denom --> variable with the current value of 1/(rij*sqrt(rij).
c It is used because the dmmsmsa subroutine is used to simplify
c the code. dmmsmsa is called for i = constant and j = 1, n
c (the case i = j is ignored because the particle cannot create
c a force on itself). Each j gives us the interaction between
c particle i and j, in the subroutine the sum of the interactions
c is place into f(1,i). In denom rij = the dot product of the
c rij vector with itself (this is known 1/r^2 term of the coulomb
c force). sqrt(rij) is the magnitude of rij. It is contained
c in denom because when we add f(1,i) to rij we want the vector
c coulomb force which is defined as 1/rij^2 * rij (unit vector).
c The negative sign in denom ensures the force vector is in the
c correct direction for a repulsive force (along rij, pointing
c away from particle j).
c
c rijmagsq --> variable containing the rij dot rij values.
c
c ddotprod --> external function that returns the dot product of
c two column vectors sent to it. (for additional info on ddotprod
c and dmmsmsa, see ddotprod.f and dmmsmsa.f)
c
c Integers -->
c
c n --> the number of particles for which coulomb forces are
c to be calculated. 2nd dimension of both f and r. This
c variable is sent to the subroutine in the subroutine call.
c
c i, j --> counter variables
c
c Declare Variables:
c----------------------------------------------------------------------
real*8 r(3,n), f(3,n)
integer n
real*8 rij(3), denom, rijmagsq, ddotprod
integer i, j
c----------------------------------------------------------------------
c intialize f elements to zero
c----------------------------------------------------------------------
do j = 1, n
do i = 1,3
f(i,j) = 0
end do
end do
c----------------------------------------------------------------------
c Main loop. Note: no calculation for i = j
c First dmmsmsa call finds ri - rj = rij. The if statement
c prevents division by zero. if ri=rj then rij = 0 then the total
c force vector will be undefined. If rij = 0 then the subroutine
c goes to 900, prints an error, and stops execution in the main
c program.
c Calculate rijmagsq and denom. The second dmmsmsa call, adds the
c current force f(1,j) to the total force f(1,i).
c----------------------------------------------------------------------
do i = 1, n
do j = 1, n
If (i .ne. j) then
call dmmsmsa(r(1,i), r(1,j), rij, 1.0d0, -1.0d0, 3, 1)
if (rij(1) .eq. 0 .and. rij(2) .eq. 0 .and. rij(3) .eq.
& 0) goto 900
rijmagsq = ddotprod(rij, rij, 3)
denom = -1.0d0 / ( rijmagsq * sqrt(rijmagsq) )
call dmmsmsa(f(1,j), rij, f(1,j), 1.0d0, denom, 3, 1)
end if
end do
end do
c----------------------------------------------------------------------
c If everything occurs normally, return to the main program
c----------------------------------------------------------------------
return
c----------------------------------------------------------------------
c Print error message for division by zero (bullet-proofing) and stop
c execution of the main program.
c----------------------------------------------------------------------
900 continue
write(0,*) 'Error (force.f): identical position vectors '
stop
end
c**********************************************************************
c**********************************************************************
c**********************************************************************
subroutine tangent(v,r,n)
implicit none
c-------------------------------------------------------------------
c Subroutine tangent:
c
c Returns the component of v tangent to r. Working with position
c vectors it returns the vector tangent to the sphere.
c
c Variables:
c
c real*8 -->
c
c v --> vector, sent by the user, that will be returned minus
c its radial component.
c
c r --> (position) vector sent by the user.
c
c proj --> takes on the values of the projection of v on r.
c This is the component we wish to remove.
c
c ddotprod --> function that returns the dot product of two
c vectors sent in the function call. (see ddotprod.f)
c
c temp --> stores the values of the r vector. Afterwards, we
c manipulate temp instead of r. This prevents the function
c from returning a different vector (r) when it returns to
c the main program.
c
c vr and rr --> take on the values of v dot r and r dot r,
c respectively. The vector tangent to r is defined as
c vtan = v - (vr/rr) * r where all variables, except those in
c parenthesis are vectors.
c
c Integers -->
c
c maxn --> the maximum number of elements in temp. We cannot
c explicitely define temp to be length n, since it was not
c given in the function call.
c
c n --> total number of elements in the v and r vectors.
c
c i --> counter variable.
c
c Vector addition is carried out using the subroutine dmmsmsa.
c For more information see dmmsmsa.f.
c
c Declare variables:
c------------------------------------------------------------------
integer maxn
parameter (maxn = 100 000)
real*8 v(n), r(n), proj, ddotprod, temp(maxn)
real*8 vr, rr
integer n, i
c------------------------------------------------------------------
c Set temp = r.
c------------------------------------------------------------------
do i = 1, n
temp(i) = r(i)
end do
c------------------------------------------------------------------
c Calculate vr, rr, and proj. notice the negative sign in proj.
c It is there so that subtraction occurs between v and temp.
c The negative sign could also appear in the dmmsmsa call.
c------------------------------------------------------------------
vr = ddotprod(v,temp,n)
rr = ddotprod(temp,temp,n)
proj = -(vr/rr)*(1.0d0)
call dmmsmsa(v,temp,v,1.0d0,proj,3,1)
c------------------------------------------------------------------
c Return to the main program.
c------------------------------------------------------------------
return
end
c**********************************************************************
c**********************************************************************
c**********************************************************************
subroutine dmmsmsa(a1,a2,a3,s1,s2,d1,d2)
implicit none
c----------------------------------------------------------------
c This subroutine returns the sum of two scaled matrices.
c
c Variables:
c
c real*8 -->
c
c a1, a2 --> matrices to be added. Both a2 and a1 are sent by
c the user.
c
c a3 --> the sum matrix returned to the main program.
c
c s1 --> user defined scalar that multiplies a1.
c
c s2 --> user defined scalar that multiplies a2.
c
c Integers -->
c
c i, j --> counter variables.
c
c d1 --) leading dimension of a1, a2, and a3. # of rows.
c
c d2 --> second dimension of a1, a2, and a3. # of columns.
c
c Declare variables:
c----------------------------------------------------------------
integer i, j, d1, d2
real*8 a1(d1,d2), a2(d1,d2), a3(d1,d2)
real*8 s1, s2
c------------------------------------------------------------------
c calculate a3 = a1*s1 + a2*s2 and return.
c------------------------------------------------------------------
do j = 1, d2
do i = 1, d1
a3(i,j) = a1(i,j) * s1 + a2(i,j) * s2
end do
end do
return
end
c**********************************************************************
c**********************************************************************
c**********************************************************************
real*8 function ddotprod(v1,v2,n)
implicit none
c------------------------------------------------------------------
c ddotprod returns the dot product of the vectors v1 and v2 sent
c by the user in the call to the function.
c
c Variables:
c
c real*8 -->
c
c v1 and v2 --> vectors that this function takes the dot product
c of. Vectors are of length n.
c
c temp --> The dot product is defined as v1x*v2x + v1y*v2y +
c v1z*v2z. temp(i) stores v1(i)*v2(i). In a loop we will
c sum the temp(i)'s into a total dotproduct.
c
c Integers -->
c
c n --> number of elements in v1, v2, and temp column vectors.
c
c i --> counter variable
c------------------------------------------------------------------
integer n, i
real*8 v1(n), v2(n), temp(n)
c------------------------------------------------------------------
c find temp. temp(i) = the product of the corresponding elements
c in v1 and v2.
c------------------------------------------------------------------
do i = 1, n
temp(i) = v1(i) * v2(i)
end do
ddotprod = 0
c------------------------------------------------------------------
c Sum the individual temp(i)'s into ddotprod and return to the
c main program.
c------------------------------------------------------------------
do i = 1, n
ddotprod = ddotprod + temp(i)
end do
c------------------------------------------------------------------
c return to the main program.
c------------------------------------------------------------------
return
end
c**********************************************************************
c**********************************************************************
c**********************************************************************
subroutine gennpart(v,n)
implicit none
c---------------------------------------------------------------
c This subroutine generates the position vectors of the n
c n particles that are to be placed on a sphere of radius r.
c It uses the rand() function to do so.
c
c Variables:
c
c real*8 -->
c
c rand --> unix function. Generates a pseudo random number
c between 0 and 1.
c
c v --> 2-d array that is returned to the main program containing
c random position vectors. Note: these position vectors are
c not normalized to the surface of the sphere.
c
c rnd --> column vector of length 3*n. It contains random values
c for each element of v. Note: it is given the value of
c rand() - rand() to give each element random sign as well as
c random value.
c
c Integers -->
c
c i, j --> counter variables.
c
c n --> second dimension of v array. Also, describes the
c total length of rnd.
c
c
c Note!!!: The case where v(1,i) = v(1,j) is possible (if not
c extremly improbable). When used with the force subroutine
c this case is screened. In other applications, the user is
c warned that this case can occur.
c
c Declare Variables:
c---------------------------------------------------------------
real*8 rand
real*8 v(3,n), rnd(3*n)
integer i, j, n, r
c---------------------------------------------------------------
c Give rnd its random values
c---------------------------------------------------------------
do i = 1 , 3*n
rnd(i) = (rand() - rand())
end do
c---------------------------------------------------------------
c Create random v array. j = 0, n-1 ; i = 1, 3, ensures that
c each value in v is unique.
c---------------------------------------------------------------
do j = 0 , n-1
do i = 1 , 3
v(i,j+1) = rnd(3*j+i)
end do
end do
c---------------------------------------------------------------
c return to the main program.
c---------------------------------------------------------------
return
end
c**********************************************************************
c**********************************************************************
c**********************************************************************
subroutine reflonsphere(v,n,r)
implicit none
c---------------------------------------------------------------
c This subroutine takes the calculated vector of an individual
c particle and changes its length to that of the radius of the
c sphere. The radius is sent in the subroutine call, in the
c main program.
c
c Variables:
c
c real*8 -->
c
c v(3,n) --> 2-d array containing n position vectors,
c corresponding to n particles. It is sent by the main
c program, the extra radial peices are removed, and it is
c sent back to the main program.
c
c ddotprod --> (function) returns the dot product of two
c column vectors sent to the function.
c
c gamma --> factor that when multiplied by the original
c vector, gives the vector magnitude r.
c
c Integers -->
c
c j --> counter variable
c
c n --> second dimension of v. Number of position vectors
c to be reflected.
c
c r --> radius of the sphere. Sent by user in subroutine
c call.
c
c Note: dmmsmsa subroutine is used to multiple v(1,i) by its
c corresponding gamma. (for more info see dmmsmsa.f)
c
c Declare variables:
c---------------------------------------------------------------
real*8 v(3,n), ddotprod
real*8 gamma
integer j, n, r
do j = 1 , n
c---------------------------------------------------------------
c sqrt(r^2/(x^2 + y^2 + z^2)) is the factor that must be
c multiplied by each component in order for the vector to have
c magnitude of r. Send v(1,j) to dmmsmsa to multiply it by
c gamma.
c---------------------------------------------------------------
gamma = sqrt(r**2 / ddotprod(v(1,j),v(1,j),3) )*1.0d0
call dmmsmsa(v(1,j),v(1,j),v(1,j),0.0d0,gamma,3,1)
end do
c---------------------------------------------------------------
c Return to main program
c---------------------------------------------------------------
return
end
c**********************************************************************
c**********************************************************************
c**********************************************************************
subroutine dumpvec(quant,npart,label,unit)
implicit none
c------------------------------------------------------------------
c Subroutine authored by Matthew W. Choptuik. When called it
c writes an entire vector to stderror. output is formatted.
c------------------------------------------------------------------
integer npart, unit
real*8 quant(3,npart)
character*(*) label
integer i, j
write(unit,1000) label
do i = 1 , npart
write(unit,1100) i, ( quant(j,i) , j = 1 , 3 )
end do
return
1000 format(/' dumpvec: ',a)
1100 format(' Particle: ',i3,1p,5e15.3,0p)
end
c**********************************************************************
c**********************************************************************
c**********************************************************************
c------------------------------------------------------------------
c C subroutine sphplot. This subroutine was authored by Dr.
c Matthew W. Choptuik. It was instrumental in the visualization
c of equilibrium configurations. (It allowed me to see the steps
c leading up to the final configuration and the final configuration
c itself)
c------------------------------------------------------------------
#include <stdio.h>
#include <gl.h>
#include <device.h>
#include <gl/sphere.h>
#include "sphplot.h"
int Ncall = 0;
int Winid;
/* If you don't like purplish spheres, try setting Sphere_color to
RED, WHITE, BLUE, GREEN, YELLOW, or CYAN */
int Sphere_color = MAGENTA;
void sphplot_(double *pos,int *pncharge,double *pradius) {
sphplot(pos,*pncharge,*pradius);
}
/* Call this routine with ncharge > 0 for normal operation,
ncharge < 0 for tumble mode until escape key is depressed
in graphics window ... */
void sphplot(double *pos,int ncharge,double radius) {
int i;
int rotx, roty, rotz;
int drotx = 2, droty = 3, drotz = 4;
long dev;
short val;
if( Ncall == 0 ) {
/* Open window, set graphics attributes and clear to black ... */
keepaspect(1,1);
foreground();
Winid = winopen("sphplot");
cmode();
doublebuffer();
gconfig();
perspective(250,1.0,0.0001,100000.0);
color(BLACK);
clear();
swapbuffers();
translate(0.0,0.0,-7.5);
/* Set sphere attributes ... */
sphmode(SPH_TESS,SPH_OCT);
sphmode(SPH_DEPTH,5);
sphmode(SPH_PRIM,SPH_POLY);
/* Queue escape-key ... */
qdevice(ESCKEY);
}
if( ncharge > 0 ) {
/* Regular draw ... */
draw_spheres(pos,ncharge,radius,0,0,0);
} else {
/* Tumble and wait for escape key ... */
rotx = 0; roty = 0; rotz = 0;
while( 1 ) {
if( qtest() ) {
dev = qread(&val);
if( (dev == ESCKEY) && val ) return;
}
draw_spheres(pos,-ncharge,radius,rotx,roty,rotz);
rotx = (rotx + drotx) % 3600;
roty = (roty + droty) % 3600;
rotz = (rotz + drotz) % 3600;
}
}
Ncall++;
}
void draw_spheres(double *pos,int ncharge,double radius,
int rotx,int roty,int rotz) {
float sphparams[4];
int i;
sphparams[3] = radius;
pushmatrix();
color(BLACK);
clear();
rotate(rotx,'x');
rotate(roty,'y');
rotate(rotz,'z');
color(Sphere_color);
for( i = 0; i < ncharge; i++ ) {
sphparams[0] = pos[3 * i];
sphparams[1] = pos[3 * i + 1];
sphparams[2] = pos[3 * i + 2];
sphdraw(sphparams);
}
swapbuffers();
popmatrix();
return;
}