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,' <>' 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,*) '<>' 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 [ '// & ']' 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 #include #include #include #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; }