/* Simple Monte Carlo simulation of a Lennard-Jones fluid Copyright 2011-2014 Cheng Zhang This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License can be found in . */ import java.applet.*; import java.awt.*; import java.awt.event.*; import java.awt.geom.*; import java.awt.image.*; import static java.lang.Math.*; import java.text.*; import javax.swing.*; import java.util.Hashtable; import java.util.Random; public class LJ3MCApp extends JApplet implements ActionListener { double rhoInit = 0.8; LJ3MC mc = new LJ3MC(rhoInit); int delay = 100; int speed = 1000; // mc steps per frame Timer timer; XYZCanvasMC canvas; // 3D drawing here JPanel cpnl, spnl; JTextField tNum = new JTextField("" + mc.N); JTextField tTemp = new JTextField("1.0"); JTextField tRho = new JTextField("" + rhoInit); JTextField tSpeed = new JTextField("" + speed); JTextField tMCAmp = new JTextField("" + mc.mcAmp); JToggleButton bStart = new JToggleButton("Start", false); JButton bReset = new JButton("Reset"); JButton bRetime = new JButton("Reset time"); JLabel lStatus = new JLabel("Status"); JTextField tAvU = new JTextField(" 0"); JTextField tAvp = new JTextField(" 0"); JTextField tAcc = new JTextField(" 0"); public static final long serialVersionUID = 1L; public void init() { tNum.setHorizontalAlignment(JTextField.CENTER); tTemp.setHorizontalAlignment(JTextField.CENTER); tRho.setHorizontalAlignment(JTextField.CENTER); tSpeed.setHorizontalAlignment(JTextField.CENTER); tMCAmp.setHorizontalAlignment(JTextField.CENTER); tAvU.setHorizontalAlignment(JTextField.RIGHT); tAvp.setHorizontalAlignment(JTextField.RIGHT); tAcc.setHorizontalAlignment(JTextField.RIGHT); Container box = getContentPane(); box.setLayout(new BorderLayout()); cpnl = new JPanel(); // create a panel for controls cpnl.setLayout(new GridLayout(19, 2)); box.add(cpnl, BorderLayout.EAST); // add controls cpnl.add(bStart); bStart.addActionListener(this); cpnl.add(bReset); bReset.addActionListener(this); cpnl.add(new JLabel(" N:")); tNum.addActionListener(this); cpnl.add(tNum); cpnl.add(new JLabel(" Temperature:")); tTemp.addActionListener(this); cpnl.add(tTemp); cpnl.add(new JLabel(" Density:")); tRho.addActionListener(this); cpnl.add(tRho); cpnl.add(new JLabel(" Steps/frame:")); tSpeed.addActionListener(this); cpnl.add(tSpeed); cpnl.add(new JLabel(" Move \u0394X:")); tMCAmp.addActionListener(this); cpnl.add(tMCAmp); cpnl.add(new JLabel(" < U/N >:")); tAvU.setEditable(false); cpnl.add(tAvU); cpnl.add(new JLabel(" < pressure >:")); tAvp.setEditable(false); cpnl.add(tAvp); cpnl.add(new JLabel(" Acc. ratio:")); tAcc.setEditable(false); cpnl.add(tAcc); cpnl.add(bRetime); bRetime.addActionListener(this); spnl = new JPanel(); // create a panel for status box.add(spnl, BorderLayout.SOUTH); lStatus.setFont(new Font("Courier", 0, 12)); spnl.add(lStatus); canvas = new XYZCanvasMC(); box.add(canvas, BorderLayout.CENTER); timer = new Timer(delay, this); timer.start(); timer.stop(); } public void start() { } public void update(Graphics g) { paint(g); } DecimalFormat df = new DecimalFormat("###0.000"); //int frame = 0; public void paint(Graphics g) { //System.out.println("frame: " + (frame++)); lStatus.setText("n = " + mc.step + ", " + "N = " + mc.N + ", " + "U/N = " + df.format(mc.U/mc.N) + ", " + "p = " + df.format(mc.p) + ";"); tAvU.setText(df.format(mc.avU.getAve()/mc.N) + " "); tAvp.setText(df.format(mc.avp.getAve()) + " "); tAcc.setText(df.format(mc.mcTot > 0 ? 1.0*mc.mcAcc/mc.mcTot : 0.0) + " "); refreshCanvas(false); cpnl.repaint(); spnl.repaint(); } public void actionPerformed(ActionEvent e) { Object src = e.getSource(); if (src == timer) { for (int i = 0; i < speed; i++) // sample a few steps mc.metropolis(); repaint(); return; } boolean adjCanvasScale = false; if (src == tTemp || src == bReset) { double kT = Double.parseDouble(tTemp.getText().trim()); if (kT < 1e-8) { kT = 1e-8; tTemp.setText("" + kT); } mc.kT = kT; mc.clearData(); } if (src == tRho || src == bReset) { double rho = Double.parseDouble( tRho.getText().trim() ); if (rho < 1e-3) { rho = 1e-3; tRho.setText("" + rho); } if (rho > 1.2) { rho = 1.2; tRho.setText("" + rho); } mc.setDensity(rho); mc.clearData(); adjCanvasScale = true; } if (src == tSpeed || src == bReset) { speed = Integer.parseInt(tSpeed.getText().trim()); if (speed < 1) { speed = 1; tSpeed.setText("" +speed); } } if (src == tMCAmp || src == bReset) { double amp = Double.parseDouble(tMCAmp.getText().trim()); if (amp < 0.) { amp = 0; tMCAmp.setText("" + amp); } if (amp > mc.L) { amp = mc.L; tMCAmp.setText("" + amp); } mc.mcAmp = amp; mc.clearData(); } if (src == bRetime) { mc.clearData(); } if (src == bStart) { boolean on = bStart.isSelected(); if (on) { timer.restart(); bStart.setText("Pause"); } else { timer.stop(); bStart.setText("Resume"); } } if (src == tNum) { int n = Integer.parseInt(tNum.getText().trim()); if (n < 2) { n = 2; tNum.setText(" " + n); } mc.N = n; mc.init(mc.rho); adjCanvasScale = true; } if (src == bReset) { if (timer.isRunning()) timer.stop(); bStart.setSelected(false); bStart.setText("Start"); mc.init(mc.rho); } refreshCanvas(adjCanvasScale); repaint(); } private void refreshCanvas(boolean adjScale) { canvas.refresh(mc.getXWrap(), mc.N, mc.xo, mc.xi, true, mc.moveAtom, mc.moveAcc, adjScale); } } class LJ3MC { static final int D = 3; public int N = 256; // number of particles 8^3/2 double rho = 0.2; double dt = 0.002; // time step for integrating Newton's equation public double L = 10.0; // side of the box double rc = 2.5; public double kT = 1.0; // temperature times Boltzmann constant double x[][]; // position, from 0 to 1 double xi[] = new double [D], xo[] = new double [D]; double xwrap[][]; // wrapped coordinates double U; // kinetic, potential, and total energy double Vir, p; // virial and pressure double Utail; // potential energy tail correction double ptail; // pressure tail correction int step; // simulation steps double mcAmp = 0.15; // amplitude of randomly moving a particle long mcAcc = 0, mcTot = 0; Ave avU = new Ave(), avp = new Ave(); /** constructor */ LJ3MC(double den) { init(den); } /** initialize system */ public void init(double den) { int i, j, k, id; step = 0; x = new double[N][D]; xwrap = new double[N][D]; int N1 = (int) (pow(2*N, 1.0/D) + .999999); // # of particles per side (cubic lattic) double a = 1./N1; for (id = 0, i = 0; i < N1 && id < N; i++) // place particles on a cubic lattice for (j = 0; j < N1 && id < N; j++) for (k = 0; k < N1 && id < N; k++) { if ((i + j + k) % 2 != 0) continue; x[id][0] = a * (i + .5); x[id][1] = a * (j + .5); x[id][2] = a * (k + .5); id++; } setDensity(den); } /** set density and recalculate tail corrections */ void setDensity(double den) { rho = den; L = pow(N/rho, 1.0/D); /* compute shifts and tail corrections */ if ((rc = 2.5) > .5*L) rc = .5*L; double irc = 1./rc, irc3 = irc*irc*irc, irc6 = irc3*irc3; Utail = 8.0*PI/9*rho*N*(irc6 - 3)*irc3; ptail = 16.0*PI/9*rho*rho*(irc6 - 1.5)*irc3; clearData(); energy(); /* re-compute energy */ } void clearData() { step = 0; mcAcc = mcTot = 0; avU.clear(); avp.clear(); } public int moveAtom = -1; public boolean moveAcc = false; Random rng = new Random(); /** Metropolis algorithm */ double xij[] = new double[D]; void metropolis() { double dU = 0, dVir = 0, r2, invr2, invr6, fscal, amp; int i, j, d; /* randomly move a particle */ moveAtom = i = (int) (rng.nextDouble() * N); amp = mcAmp / L; for (d = 0; d < D; d++) { xo[d] = x[i][d]; xi[d] = xo[d] + amp * (rng.nextDouble()*2.0 - 1); } for (j = 0; j < N; j++) { if (i == j) continue; /* old energy */ rv3Diff(xij, x[i], x[j]); vpbc(xij); r2 = vsqr(xij); if (r2 < rc*rc) { invr2 = 1.0/r2; invr6 = invr2 * invr2 * invr2; fscal = invr6 * (48.0 * invr6 - 24.0); dVir -= fscal; dU -= 4.0 * invr6 * (invr6 - 1.0); } /* new energy */ rv3Diff(xij, xi, x[j]); vpbc(xij); r2 = vsqr(xij); if (r2 < rc*rc) { invr2 = 1.0/r2; invr6 = invr2 * invr2 * invr2; fscal = invr6 * (48.0 * invr6 - 24.0); dVir += fscal; dU += 4.0 * invr6 * (invr6 - 1.0); } } mcTot++; if (dU <= 0 || rng.nextDouble() < exp(-dU/kT) ) { mcAcc++; for (d = 0; d < D; d++) x[i][d] = xi[d]; U += dU; Vir += dVir; calcPressure(); moveAcc = true; //System.out.println("a U " + U + " dU " + dU); //energy(); //System.out.println("b U " + U); } else { moveAcc = false; } step++; avU.add(U); avp.add(p); } /** Remove the center of mass motion */ void rmcom() { for (int k = 0; k < D; k++) { double xc = 0; for (int i = 0; i < N; i++) xc += x[i][k]; xc /= N; for (int i = 0; i < N; i++) x[i][k] -= xc; xi[k] -= xc; xo[k] -= xc; } } /** calculate potential energy and force * half-box cutoff, minimal distance image */ void energy() { int i, j, d; double invr2, invr6, r2, fscal; // clear the energy and pressure U = Vir = 0.0; // loop over pairs of particles and compute energy for (i = 0; i < N - 1; i++) { for (j = i + 1; j < N; j++) { rv3Diff(xij, x[i], x[j]); vpbc(xij); r2 = vsqr(xij); if (r2 > rc*rc) continue; invr2 = 1.0/r2; invr6 = invr2 * invr2 * invr2; fscal = invr6 * (48.0 * invr6 - 24.0); Vir += fscal; U += 4.0 * invr6 * (invr6 - 1.0); } } U += Utail; calcPressure(); } void calcPressure() { p = rho*kT + Vir/(3*L*L*L) + ptail; } private void rv3Diff(double c[], double a[], double b[]) { c[0] = a[0] - b[0]; c[1] = a[1] - b[1]; c[2] = a[2] - b[2]; } /** image with the minimal distance */ double pbc(double a) { return L * (a - ((int)(a + 1000.5) - 1000)); } void vpbc(double v[]) { v[0] = pbc(v[0]); v[1] = pbc(v[1]); v[2] = pbc(v[2]); } double vsqr(double v[]) { return v[0]*v[0] + v[1]*v[1] + v[2]*v[2]; } /** wrap coordinates back into box */ public double [][] getXWrap() { rmcom(); double xx; for (int i = 0; i < N; i++) for (int d = 0; d < D; d++) { xx = x[i][d] + 1000.0; xwrap[i][d] = (xx - (int) xx) * L; } for (int d = 0; d < D; d++) { xx = xi[d] + 1000.0; xi[d] = (xx - (int) xx) * L; xx = xo[d] + 1000.0; xo[d] = (xx - (int) xx) * L; } return xwrap; } } /** Average and standard deviation */ class Ave { double cnt, xsm, x2sm; public void clear() { cnt = xsm = x2sm = 0; } public void add(double x) { cnt += 1; xsm += x; x2sm += x*x; } public double getAve() { return (cnt > 0) ? xsm/cnt : 0; } public double getVar() { if (cnt <= 1) return 0; double av = xsm/cnt; return x2sm/cnt - av*av; } public double getStd() { return Math.sqrt( getVar() ); } } /** Panel for drawing particles */ class XYZCanvas extends JPanel implements MouseListener, MouseMotionListener, MouseWheelListener { Image img; Graphics imgG; Dimension imgSize; double realSpan; // size of a real span (saved as a copy) double zoomScale = 1.0; public Matrix3D viewMatrix = new Matrix3D(); // view matrix private Matrix3D tmpMatrix = new Matrix3D(); // temporary matrix int mouseX, mouseY; // mouse position XYZModel model; public static final long serialVersionUID = 2L; public XYZCanvas() { super(); addMouseListener(this); addMouseMotionListener(this); addMouseWheelListener(this); } /** Prepare a buffer for the image, return if the canvas is ready * Only need to call this when the size of the canvas is changed, * Since we automatically detect the size change in paint() * only call this on start */ public boolean newImgBuf() { Dimension sz = getSize(); if (sz.width == 0 || sz.height == 0) return false; // quit if the current image already has the right size if (img != null && imgG != null && sz.equals(imgSize)) return true; img = createImage(sz.width, sz.height); if (imgG != null) imgG.dispose(); imgG = img.getGraphics(); imgSize = sz; return true; } public void update(Graphics g) { if (img == null) g.clearRect(0, 0, getSize().width, getSize().height); paintComponent(g); } protected void paintComponent(Graphics g) { super.paintComponent(g); if (model != null) { newImgBuf(); // refresh the image buffer if necessary // compute the real-to-screen ratio, this variable differs // from model.real2Screen by zoomScale Dimension sz = getSize(); double real2Screen0 = model.getScaleFromSpan(realSpan, sz.width, sz.height); model.setMatrix(viewMatrix, real2Screen0 * zoomScale, sz.width/2, sz.height/2); imgG.setColor(Color.BLACK); imgG.fillRect(0, 0, sz.width, sz.height); model.paint(imgG); g.drawImage(img, 0, 0, this); } } /** Refresh the coordinates * x[][] is the wrapped coordinates * `n' may be less than x.length */ public void refresh(double x[][], int n, boolean center, boolean adjScale) { if (model == null) { model = new XYZModel(); adjScale = true; } model.updateXYZ(x, n, center); if ( adjScale ) realSpan = model.getSpan(x, n); repaint(); } /** Event handling */ public void mouseClicked(MouseEvent e) { } public void mousePressed(MouseEvent e) { mouseX = e.getX(); mouseY = e.getY(); // consume this event so that it will not be processed // in the default manner e.consume(); } public void mouseReleased(MouseEvent e) { } public void mouseEntered(MouseEvent e) { } public void mouseExited(MouseEvent e) { } public void mouseDragged(MouseEvent e) { int x = e.getX(); int y = e.getY(); tmpMatrix.unit(); tmpMatrix.xrot(360.0 * (mouseY - y) / getSize().height); tmpMatrix.yrot(360.0 * (x - mouseX) / getSize().width); viewMatrix.mult(tmpMatrix); repaint(); mouseX = x; mouseY = y; e.consume(); } public void mouseMoved(MouseEvent e) { } public void mouseWheelMoved(MouseWheelEvent e) { int notches = e.getWheelRotation(); if ((zoomScale -= 0.05f*notches) < 0.09999f) zoomScale = 0.1f; repaint(); } } /** Panel for drawing particles */ class XYZCanvasMC extends XYZCanvas { public static final long serialVersionUID = 3L; /** Get 3D position * `n' may be less than x.length */ private double [][] appendXMove(double x[][], int n, double xMove[]) { int D = x[0].length; double r[][] = new double[n + 1][3]; // copy the direct coordinates for (int i = 0; i < n; i++) { r[i][0] = x[i][0]; r[i][1] = x[i][1]; r[i][2] = (D == 2) ? 0.0 : x[i][2]; } if (xMove == null) xMove = x[0]; // the last is the trial position r[n][0] = xMove[0]; r[n][1] = xMove[1]; r[n][2] = (D == 2) ? 0 : xMove[2]; return r; } /** Refresh the coordinates * x[][] is the wrapped coordinates * `n' may be less than x.length */ public void refresh(double x[][], int n, double xo[], double xi[], boolean center, int mvAtom, boolean mvAcc, boolean adjScale) { if (model == null) { model = new XYZModelMC(); adjScale = true; } XYZModelMC modelMC = (XYZModelMC) model; double xyz[][] = appendXMove(x, n, mvAcc ? xo : xi); modelMC.updateXYZ(xyz, n + 1, center); modelMC.setMoveAtom(mvAtom, mvAcc); if ( adjScale ) realSpan = model.getSpan(x, n); repaint(); } } /** A set of atoms with 3D coordinates */ class XYZModel { double realXYZ[][]; // 3D real coordinates [np][3] int screenXYZ[][]; // 3D screen coordinates in pixels [np][3] // only the first two dimensions are used in drawing int zOrder[]; // z-order, larger is closer to the viewer int np = -1; boolean transformed; // rotation/scaling/translation matrix for the conversion from realXYZ[] to screenXYZ[] Matrix3D mat = new Matrix3D(); double real2Screen = 50.0; // real size --> screen size double ballSize = 0.5; // ball size (radius) in terms of the real coordinates // 0.5 for hard spheres Atom atomDefault = new Atom(0.1, 0.7, 0.1, 1.0, 0.5, 0.5, 1.0); Atom atoms[]; // for colors of atoms XYZModel() {} /** Set the color of a particular atom */ void setAtom(int i, Atom atom) { if ( i >= 0 && i < atoms.length ) atoms[i] = atom; } /** Refresh coordinates * x[0..n-1][3] is a three-dimensional vector * n can be less than x.length */ void updateXYZ(double x[][], int n, boolean center) { if (n != np) { //System.out.printf("XYZModel.updateXYZ n %d --> %d, (%g, %g, %g) (%g, %g, %g)\n", np, n, x[0][0], x[0][1], x[0][2], x[n-1][0], x[n-1][1], x[n-1][2]); np = n; realXYZ = new double [np][3]; screenXYZ = new int [np][3]; zOrder = new int [np]; atoms = new Atom [np]; for (int i = 0; i < np; i++) atoms[i] = atomDefault; } for (int d = 0; d < 3; d++) { double xc = 0; if ( center ) { for (int i = 0; i < np; i++) xc += x[i][d]; xc /= np; } for (int i = 0; i < np; i++) realXYZ[i][d] = x[i][d] - xc; } transformed = false; } /** Set the view matrix * `s' is the scaling factor of translating real coordinates * to the screen coordinates * (x0, y0) the screen coordinates of the center */ void setMatrix(Matrix3D viewMat, double s, double x0, double y0) { mat.unit(); mat.mult(viewMat); mat.scale(s, s, s); real2Screen = s; mat.translate(x0, y0, 0); transformed = false; } /** Get the span of the model * `n' may be less than x.length */ double getSpan(double x[][], int n) { int dim = x[0].length; double realSpan = 0, del, fw, fh; for (int d = 0; d < dim; d++) { double xmin = 1e30, xmax = -1e30; for (int i = 0; i < n; i++) if ( x[i][d] < xmin ) xmin = x[i][d]; else if ( x[i][d] > xmax ) xmax = x[i][d]; if ( (del = xmax - xmin) > realSpan ) realSpan = del; } return realSpan; } /** Translate a given span, return the real-to-screen ratio * `w' and `h' are the width and height of the screen in pixels */ double getScaleFromSpan(double realSpan, int w, int h) { realSpan += ballSize * 2; // add two radii double fw = w / realSpan; double fh = h / realSpan; double facShrink = 0.9; // shrink a bit for the margin return (fw < fh ? fw : fh) * facShrink; } /** Compute the Z-order */ void getZOrder() { // transform the coordinates if ( !transformed ) { mat.transform(realXYZ, screenXYZ, np); transformed = true; } // bubble sort z-order // zOrder[0] is the fartherest from the viewer // zOrder[np - 1] is the nearest to the viewer for (int i = 0; i < np; i++) zOrder[i] = i; for (int i = 0; i < np; i++) { // find the particle with the smallest z int jm = i, k; for (int j = i + 1; j < np; j++) if (screenXYZ[zOrder[j]][2] < screenXYZ[zOrder[jm]][2]) jm = j; if (jm != i) { k = zOrder[i]; zOrder[i] = zOrder[jm]; zOrder[jm] = k; } } } /** Draw atom */ void drawAtom(Graphics g, int iz) { int i = zOrder[iz]; Atom atom = atoms[i]; if (atom == null) return; double zMin = screenXYZ[zOrder[0]][2]; // fartherest from the viewer double zMax = screenXYZ[zOrder[np - 1]][2]; // closest to the viewer int greyScale = Atom.nBalls - 1; if (zMin != zMax) // zMin == zMax means the two-dimensional case greyScale = (int) (Atom.nBalls * (screenXYZ[i][2] - zMin) / (zMax - zMin) - 1e-6); // the atom closest to the viewer has a greyScale of Atom.nBalls - 1 // the atom fartherest from the viewer has a greyScale of 0 double radius = ballSize * atom.relRadius * real2Screen; atom.paint(g, screenXYZ[i][0], screenXYZ[i][1], greyScale, radius); } /** Paint this model to the graphics context `g' */ void paint(Graphics g) { if (realXYZ == null || np <= 0) return; getZOrder(); for (int iz = 0; iz < np; iz++) drawAtom(g, iz); } } /** XYZModel for a typical MC simulation */ class XYZModelMC extends XYZModel { int moveAtom = -1; boolean moveAcc = false; // failed move Atom atomTrial = new Atom(1.0, 0.1, 0.1, 1.0, 0.5, 0.3, 2.0); // failed trial postion Atom atomFailed = new Atom(0.0, 0.0, 0.4); // new and old position // successful move Atom atomMoved = new Atom(0.1, 1.0, 0.5); // trial and new position Atom atomOldPos = new Atom(0.3, 0.3, 0.3, 1.0, 0.5, 0.5, 2.0); // old position /** Constructor from the real coordinates * x[][3] is a three-dimensional vector */ XYZModelMC() { atomDefault = new Atom(0.1, 0.2, 1.0, 1.0, 0.6, 0.4, 1.0); } /** Set the atom */ public void setMoveAtom(int mvAtom, boolean mvAcc) { moveAtom = mvAtom; moveAcc = mvAcc; if (mvAtom < 0) return; for (int i = 0; i < np; i++) atoms[i] = atomDefault; if ( moveAtom < 0 || moveAtom >= np - 1 ) { atoms[np - 1] = null; } else { atoms[moveAtom] = moveAcc ? atomMoved : atomFailed; atoms[np - 1] = moveAcc ? atomOldPos : atomTrial; } } } /** Atom: a ball */ class Atom { private final static int R = 120; private final static int hx = 45; // (hx, hy) is the offset of the spot light from the center private final static int hy = 45; private static int maxr; // maximal distance from the spotlight public final static int nBalls = 16; // shades of grey // spotlight brightness (0.0, 1.0) // 1.0 means the spotlight is pure white // 0.0 means the spotlight is the same as the solid color double spotlightAmp = .4; // color damping along the distance from the spotlight double rContrast = 0.7; // z-depth contrast (0.0, inf) // inf means the fartherest atom is completely dark // 0.0 means the fartherest atom is the same as the // nearest atom double zContrast = 2.0; private static byte data[]; static { // data[] is a bitmap image of the ball of radius R data = new byte[R * 2 * R * 2]; for (int Y = -R; Y < R; Y++) { int x0 = (int) (Math.sqrt(R * R - Y * Y) + 0.5); for (int X = -x0; X < x0; X++) { // sqrt(x^2 + y^2) gives distance from the spot light int x = X + hx, y = Y + hy; int r = (int) (Math.sqrt(x * x + y * y) + 0.5); // set the maximal intensity to the maximal distance // (in pixels) from the spot light if (r > maxr) maxr = r; data[(Y + R) * (R * 2) + (X + R)] = (r <= 0) ? 1 : (byte) r; } } } // the following variables are atom dependent private int Rl = 100, Gl = 100, Bl = 100; private Image balls[]; // 0..nBalls-1, at different z distances /** Constructor */ Atom(int r, int g, int b) { setRGB(r, g, b); } /** colors that range from 0 to 1 */ Atom(double r, double g, double b) { setRGB(r, g, b); } double relRadius = 1; // only used to save the radius Atom(double r, double g, double b, double rad) { this(r, g, b); relRadius = rad; } Atom(double r, double g, double b, double rad, double rcontrast, double spotlight, double zcontrast) { relRadius = rad; rContrast = rcontrast; spotlightAmp = spotlight; zContrast = zcontrast; setRGB(r, g, b); } /** Set color */ void setRGB(int r, int g, int b) { Rl = r; Gl = g; Bl = b; makeBalls(); } void setRGB(double r, double g, double b) { Rl = (int) (256*r - 1e-6); Gl = (int) (256*g - 1e-6); Bl = (int) (256*b - 1e-6); makeBalls(); } /** Linearly interpolate colors */ private int blend(double fg, double bg, double fgfactor) { return (int) (bg + (fg - bg) * fgfactor); } // need a component instance to call createImage() private static Component component = new Applet(); /** Prepare ball images with different sizes */ private void makeBalls() { balls = new Image[nBalls]; byte red[] = new byte[256]; byte green[] = new byte[256]; byte blue[] = new byte[256]; for (int id = 0; id < nBalls; id++) { // smaller `b' means closer to black // if id == 0 (fartherest from the viewer) // b = 1/(1 + zContrast) // the outer blend() gives a color close to bgGrey // if id == nBalls - 1 (closest to the viewer), // b = 1, the outer blend() gives the color of // the inner blend() double b = (zContrast*id/(nBalls - 1) + 1) / (zContrast + 1); for (int i = maxr; i >= 0; --i) { // closeness to the spotlight double q = 1 - 1. * i / maxr; // dampness of the color along the radius double p = 1 - rContrast * i / maxr; // contrast of the spotlight // if i == 0 (closest to the spotlight), // d = 1.0 - spotlightAmp, the inner // blend() gives a color close to 255 // (if spotlightAmp == 1). // if i == maxr (fartherest from the spotlight), // d = 1.0, the inner blend() gives // the foreground color, i.e., Rl, Gl, Bl // Thus, the inner blend() depends on the distance // from the spotlight, i == 0 means to be closest // to the spotlight double d = 1 - q * spotlightAmp; red[i] = (byte) blend(blend(Rl*p, 255, d), 0, b); green[i] = (byte) blend(blend(Gl*p, 255, d), 0, b); blue[i] = (byte) blend(blend(Bl*p, 255, d), 0, b); } // 256 color model IndexColorModel model = new IndexColorModel( 8, maxr + 1, red, green, blue, 0); balls[id] = component.createImage( new MemoryImageSource(R * 2, R * 2, model, data, 0, R * 2) ); } } /** Draw a ball at screen coordinate (x, y) with a ball index `id' * (0, 0) represents the top-left corner * `x', `y', `radius' are given in pixels * the ball index (gray code) `id' can be 0 to 15 */ void paint(Graphics gc, int x, int y, int id, double radius) { if (balls == null) makeBalls(); Image img = balls[id]; // id = [0..15] int size = (int) (radius * 2 + .5); gc.drawImage(img, x - size/2, y - size/2, size, size, null); //System.out.println("" + x + " " + y + " " + id + " " + radius); } } class Matrix3D { double xx, xy, xz, xo; double yx, yy, yz, yo; double zx, zy, zz, zo; static final double pi = 3.14159265; /** Create a new unit matrix */ Matrix3D() { xx = 1.0; yy = 1.0; zz = 1.0; } /** Scale along each axis independently */ void scale(double xf, double yf, double zf) { xx *= xf; xy *= xf; xz *= xf; xo *= xf; yx *= yf; yy *= yf; yz *= yf; yo *= yf; zx *= zf; zy *= zf; zz *= zf; zo *= zf; } /** Translate the origin */ void translate(double x, double y, double z) { xo += x; yo += y; zo += z; } /** Rotate theta degrees around the y axis */ void yrot(double theta) { theta *= (pi / 180); double ct = Math.cos(theta); double st = Math.sin(theta); double Nxx = (xx * ct + zx * st); double Nxy = (xy * ct + zy * st); double Nxz = (xz * ct + zz * st); double Nxo = (xo * ct + zo * st); double Nzx = (zx * ct - xx * st); double Nzy = (zy * ct - xy * st); double Nzz = (zz * ct - xz * st); double Nzo = (zo * ct - xo * st); xo = Nxo; xx = Nxx; xy = Nxy; xz = Nxz; zo = Nzo; zx = Nzx; zy = Nzy; zz = Nzz; } /** Rotate theta degrees about the x axis */ void xrot(double theta) { theta *= (pi / 180); double ct = Math.cos(theta); double st = Math.sin(theta); double Nyx = (yx * ct + zx * st); double Nyy = (yy * ct + zy * st); double Nyz = (yz * ct + zz * st); double Nyo = (yo * ct + zo * st); double Nzx = (zx * ct - yx * st); double Nzy = (zy * ct - yy * st); double Nzz = (zz * ct - yz * st); double Nzo = (zo * ct - yo * st); yo = Nyo; yx = Nyx; yy = Nyy; yz = Nyz; zo = Nzo; zx = Nzx; zy = Nzy; zz = Nzz; } /** Rotate theta degrees about the z axis */ void zrot(double theta) { theta *= pi / 180; double ct = Math.cos(theta); double st = Math.sin(theta); double Nyx = (yx * ct + xx * st); double Nyy = (yy * ct + xy * st); double Nyz = (yz * ct + xz * st); double Nyo = (yo * ct + xo * st); double Nxx = (xx * ct - yx * st); double Nxy = (xy * ct - yy * st); double Nxz = (xz * ct - yz * st); double Nxo = (xo * ct - yo * st); yo = Nyo; yx = Nyx; yy = Nyy; yz = Nyz; xo = Nxo; xx = Nxx; xy = Nxy; xz = Nxz; } /** Multiply this matrix by a second: M = M*R */ void mult(Matrix3D rhs) { double lxx = xx * rhs.xx + yx * rhs.xy + zx * rhs.xz; double lxy = xy * rhs.xx + yy * rhs.xy + zy * rhs.xz; double lxz = xz * rhs.xx + yz * rhs.xy + zz * rhs.xz; double lxo = xo * rhs.xx + yo * rhs.xy + zo * rhs.xz + rhs.xo; double lyx = xx * rhs.yx + yx * rhs.yy + zx * rhs.yz; double lyy = xy * rhs.yx + yy * rhs.yy + zy * rhs.yz; double lyz = xz * rhs.yx + yz * rhs.yy + zz * rhs.yz; double lyo = xo * rhs.yx + yo * rhs.yy + zo * rhs.yz + rhs.yo; double lzx = xx * rhs.zx + yx * rhs.zy + zx * rhs.zz; double lzy = xy * rhs.zx + yy * rhs.zy + zy * rhs.zz; double lzz = xz * rhs.zx + yz * rhs.zy + zz * rhs.zz; double lzo = xo * rhs.zx + yo * rhs.zy + zo * rhs.zz + rhs.zo; xx = lxx; xy = lxy; xz = lxz; xo = lxo; yx = lyx; yy = lyy; yz = lyz; yo = lyo; zx = lzx; zy = lzy; zz = lzz; zo = lzo; } /** Recover the unit matrix */ void unit() { xo = 0; xx = 1; xy = 0; xz = 0; yo = 0; yx = 0; yy = 1; yz = 0; zo = 0; zx = 0; zy = 0; zz = 1; } /** Transform np points from v into tv. * v contains the input coordinates in floating point. * Three successive entries in the array constitute a point. * tv ends up holding the transformed points as integers; * three successive entries per point */ void transform(double v[][], int tv[][], int np) { // np can be different from v.length for (int i = 0; i < np; i++) { double x = v[i][0], y = v[i][1], z = v[i][2]; tv[i][0] = (int) (xx * x + xy * y + xz * z + xo); tv[i][1] = (int) (yx * x + yy * y + yz * z + yo); tv[i][2] = (int) (zx * x + zy * y + zz * z + zo); } } }