// base.java (C) 2001 by Paul Falstad, www.falstad.com // this file must be run through cpp with one of BUILD_E, BUILD_V, // BUILD_M defined before compiling with java compiler. This was done // to reduce class file sizes. import java.io.InputStream; import java.awt.*; import java.awt.image.*; import java.applet.Applet; import java.applet.AudioClip; import java.util.Vector; import java.util.Hashtable; import java.util.Enumeration; import java.io.File; import java.net.URL; import java.util.Random; import java.awt.image.MemoryImageSource; import java.lang.Math; import java.text.NumberFormat; import java.lang.reflect.Constructor; import java.lang.reflect.Method; import java.awt.event.*; #ifdef BUILD_E #define BUILD_CASE_EMV(E,M,V) E #elif defined(BUILD_M) #define BUILD_CASE_EMV(E,M,V) M #else #define BUILD_CASE_EMV(E,M,V) V #endif class VecDemoCanvas extends Canvas { VecDemoFrame pg; VecDemoCanvas(VecDemoFrame p) { pg = p; } public Dimension getPreferredSize() { return new Dimension(300,400); } public void update(Graphics g) { pg.updateVecDemo(g); } public void paint(Graphics g) { pg.updateVecDemo(g); } }; class VecDemoLayout implements LayoutManager { public VecDemoLayout() {} public void addLayoutComponent(String name, Component c) {} public void removeLayoutComponent(Component c) {} public Dimension preferredLayoutSize(Container target) { return new Dimension(500, 500); } public Dimension minimumLayoutSize(Container target) { return new Dimension(100,100); } public void layoutContainer(Container target) { int barwidth = 0; int i; for (i = 1; i < target.getComponentCount(); i++) { Component m = target.getComponent(i); if (m.isVisible()) { Dimension d = m.getPreferredSize(); if (d.width > barwidth) barwidth = d.width; } } Insets insets = target.insets(); int targetw = target.size().width - insets.left - insets.right; int cw = targetw-barwidth; int targeth = target.size().height - (insets.top+insets.bottom); target.getComponent(0).move(insets.left, insets.top); target.getComponent(0).resize(cw, targeth); cw += insets.left; int h = insets.top; for (i = 1; i < target.getComponentCount(); i++) { Component m = target.getComponent(i); if (m.isVisible()) { Dimension d = m.getPreferredSize(); if (m instanceof DecentScrollbar || m instanceof TextComponent) d.width = barwidth; if (m instanceof Choice && d.width > barwidth) d.width = barwidth; if (m instanceof Label) { h += d.height/5; d.width = barwidth; } m.move(cw, h); m.resize(d.width, d.height); h += d.height; } } } }; public class VecDemo extends Applet implements ComponentListener { static VecDemoFrame ogf; void destroyFrame() { if (ogf != null) ogf.dispose(); ogf = null; repaint(); } boolean started = false; public void init() { addComponentListener(this); } public static void main(String args[]) { ogf = new VecDemoFrame(null); ogf.init(); } void showFrame() { if (ogf == null) { started = true; ogf = new VecDemoFrame(this); ogf.init(); repaint(); } } public void paint(Graphics g) { String s = "Applet is open in a separate window."; if (!started) s = "Applet is starting."; else if (ogf == null) s = "Applet is finished."; else ogf.show(); g.drawString(s, 10, 30); } public void componentHidden(ComponentEvent e){} public void componentMoved(ComponentEvent e){} public void componentShown(ComponentEvent e) { showFrame(); } public void componentResized(ComponentEvent e) {} public void destroy() { if (ogf != null) ogf.dispose(); ogf = null; repaint(); } }; class VecDemoFrame extends Frame implements ComponentListener, ActionListener, MouseMotionListener, MouseListener, ItemListener, DecentScrollbarListener { Thread engine = null; Dimension winSize; Rectangle viewMain, viewAxes; Image dbimage; Image backimage; MemoryImageSource imageSource; int pixels[]; VecDemo applet; Random random; public String getAppletInfo() { return "VecDemo by Paul Falstad"; } static final double pi = 3.14159265358979323846; int getrand(int x) { int q = random.nextInt(); if (q < 0) q = -q; return q % x; } VecDemoCanvas cv; Checkbox stoppedCheck; Button resetButton; Button kickButton; Checkbox reverseCheck; Button infoButton; Choice functionChooser; Choice dispChooser; static final int DISP_PART_VELOC = 0; static final int DISP_PART_FORCE = 1; static final int DISP_VECTORS = 2; static final int DISP_NONE = 3; static final int DISP_CURLERS = 4; Label partCountLabel; Label textFieldLabel; Label strengthLabel; DecentScrollbar partCountBar; DecentScrollbar strengthBar; DecentScrollbar aux1Bar; DecentScrollbar aux2Bar; DecentScrollbar aux3Bar; double fieldStrength, barFieldStrength; Color darkYellow = new Color(144, 144, 0); final double lineWidth = .001; class AuxBar { DecentScrollbar bar; Label label; AuxBar(Label l, DecentScrollbar b) { label = l; bar = b; } }; AuxBar auxBars[]; Label vecDensityLabel; DecentScrollbar vecDensityBar; Label potentialLabel; DecentScrollbar potentialBar; Label lineDensityLabel; DecentScrollbar lineDensityBar; Choice modeChooser; Choice floorColorChooser; Choice floorLineChooser; TextField textFields[]; int reverse; int xpoints[]; int ypoints[]; GridElement grid[][]; Particle particles[]; FieldVector vectors[]; int vecCount; int density[][]; Checkbox flatCheck; boolean isFlat; double viewAngle, viewAngleDragStart; double viewZoom = 1.6, viewZoomDragStart; double viewAngleCos = 1, viewAngleSin = 0; double viewHeight = 2, viewHeightDragStart; double viewDistance = 5; int integralX = -1, integralY; int vectorSpacing = 16; int currentStep; boolean showA; boolean parseError; Color fieldColors[]; static final int gridsize = 80; static final int densitygridsize = 16; static final double densitygroupsize = 2./densitygridsize; static final int maxParticleCount = 2500; boolean functionChanged; boolean backgroundChanged; boolean dragging, draggingView; int oldDragX, oldDragY, dragX, dragY, dragStartX, dragStartY; double dragZoomStart; Vector functionList; VecFunction curfunc; int pause = 20; static final int MOT_VELOCITY = 0; static final int MOT_FORCE = 1; static final int MOT_CURLERS = 2; static final int MOT_EQUIPOTENTIAL = 3; static final int FC_FIELD = 0; static final int FC_POTENTIAL = 1; static final int FC_NONE = 2; static final int FC_DIV = 3; static final int FC_CURL = 4; static final int FL_NONE = 0; static final int FL_GRID = 1; static final int FL_EQUIP = 2; static final int FL_LINES = 3; static final int MODE_VIEW_ROTATE = 0; static final int MODE_VIEW_ZOOM = 1; static final int MODE_LINE_INT = 2; static final int MODE_SURF_INT = 3; boolean useBufferedImage = false; VecDemoFrame(VecDemo a) { super(BUILD_CASE_EMV("2-D Electrostatic Fields Applet v1.4b", null, "2-D Vector Fields Applet v1.4b")); applet = a; } public void init() { try { String param = applet.getParameter("PAUSE"); if (param != null) pause = Integer.parseInt(param); } catch (Exception e) { } String jv = System.getProperty("java.class.version"); double jvf = new Double(jv).doubleValue(); if (jvf >= 48) useBufferedImage = true; functionList = new Vector(); VecFunction vf = new InverseRadial(); int ct = 0; while (vf != null) { functionList.addElement(vf); vf = vf.createNext(); if (ct == 1000) { System.out.print("setup loop\n"); return; } } Color particleColors[] = new Color[27]; int i; for (i = 0; i != 27; i++) particleColors[i] = new Color(((i%3)+1)*85, ((i/3)%3+1)*85, ((i/9)%3+1)*85); random = new Random(); particles = new Particle[maxParticleCount]; for (i = 0; i != maxParticleCount; i++) { particles[i] = new Particle(); particles[i].color = particleColors[i % 27]; } xpoints = new int[4]; ypoints = new int[4]; density = new int[densitygridsize][densitygridsize]; setLayout(new VecDemoLayout()); cv = new VecDemoCanvas(this); cv.addComponentListener(this); cv.addMouseMotionListener(this); cv.addMouseListener(this); add(cv); /*infoButton = new Button("Function Info"); add(infoButton); infoButton.addActionListener(this);*/ functionChooser = new Choice(); for (i = 0; i != functionList.size(); i++) functionChooser.add("Setup: " + ((VecFunction) functionList.elementAt(i)).getName()); add(functionChooser); functionChooser.addItemListener(this); floorColorChooser = new Choice(); floorColorChooser.add("Color: field magnitude"); floorColorChooser.add("Color: potential"); floorColorChooser.add("Color: none"); #ifdef BUILD_E floorColorChooser.add("Color: charge"); #else floorColorChooser.add("Color: divergence"); floorColorChooser.add("Color: curl z"); #endif floorColorChooser.addItemListener(this); add(floorColorChooser); floorLineChooser = new Choice(); floorLineChooser.add("Floor: no lines"); floorLineChooser.add("Floor: grid"); floorLineChooser.add("Floor: equipotentials"); #ifdef BUILD_V floorLineChooser.add("Floor: streamlines"); #else floorLineChooser.add("Floor: field lines"); #endif floorLineChooser.addItemListener(this); add(floorLineChooser); floorLineChooser.select(FL_EQUIP); flatCheck = new Checkbox("Flat View"); flatCheck.addItemListener(this); add(flatCheck); dispChooser = new Choice(); dispChooser.addItemListener(this); setupDispChooser(true); add(dispChooser); modeChooser = new Choice(); modeChooser.add("Mouse = Adjust Angle"); modeChooser.add("Mouse = Adjust Zoom"); modeChooser.add("Mouse = Line Integral"); modeChooser.add("Mouse = Surface Integral"); modeChooser.addItemListener(this); add(modeChooser); stoppedCheck = new Checkbox("Stopped"); stoppedCheck.addItemListener(this); add(stoppedCheck); reverseCheck = new Checkbox("Reverse"); reverseCheck.addItemListener(this); add(reverseCheck); resetButton = new Button("Reset"); add(resetButton); resetButton.addActionListener(this); kickButton = new Button("Kick"); add(kickButton); kickButton.addActionListener(this); kickButton.disable(); add(strengthLabel = new Label("Field Strength", Label.CENTER)); add(strengthBar = new DecentScrollbar(this, 80, 1, 120)); add(partCountLabel = new Label("Number of Particles", Label.CENTER)); add(partCountBar = new DecentScrollbar(this, 500, 1, maxParticleCount)); add(vecDensityLabel = new Label("Vector Density", Label.CENTER)); add(vecDensityBar = new DecentScrollbar(this, 32, 2, 64)); /*add(lineDensityLabel = new Label("Line Density", Label.CENTER)); add(lineDensityBar = new DecentScrollbar(this, 5, 3, 8));*/ add(potentialLabel = new Label("Potential", Label.CENTER)); add(potentialBar = new DecentScrollbar(this, 250, 0, 1000)); Label lb; auxBars = new AuxBar[3]; add(lb = new Label("Aux 1", Label.CENTER)); add(aux1Bar = new DecentScrollbar(this, 0, 0, 100)); auxBars[0] = new AuxBar(lb, aux1Bar); add(lb = new Label("Aux 2", Label.CENTER)); add(aux2Bar = new DecentScrollbar(this, 0, 0, 100)); auxBars[1] = new AuxBar(lb, aux2Bar); add(lb = new Label("Aux 3", Label.CENTER)); add(aux3Bar = new DecentScrollbar(this, 0, 0, 100)); auxBars[2] = new AuxBar(lb, aux3Bar); #ifdef BUILD_V add(textFieldLabel = new Label("", Label.CENTER)); #endif textFields = new TextField[2]; for (i = 0; i != 2; i++) { add(textFields[i] = new TextField()); textFields[i].addActionListener(this); } fieldColors = new Color[513]; int grayLevel = 76; for (i = 0; i != 256; i++) { int rb = grayLevel+(128-grayLevel)*i/255; int g = grayLevel+(255-grayLevel)*i/255; int col = (255<<24) | (g<<8) | (rb<<16) | (rb); fieldColors[i] = new Color(col); } for (i = 0; i != 256; i++) { int col = (255<<24) | (255<<8) | ((i/2+128) * (0x10001)); fieldColors[i+256] = new Color(col); } fieldColors[512] = fieldColors[511]; add(new Label("", Label.CENTER)); reinit(); cv.setBackground(Color.black); cv.setForeground(Color.lightGray); resize(650, 500); handleResize(); Dimension screen = getToolkit().getScreenSize(); Dimension x = getSize(); setLocation((screen.width - x.width)/2, (screen.height - x.height)/2); functionChanged(); dispChooserChanged(); show(); } void handleResize() { Dimension d = winSize = cv.getSize(); if (winSize.width == 0) return; dbimage = createImage(d.width, d.height); scaleworld(); viewMain = new Rectangle(winSize); viewAxes = new Rectangle(winSize.width-100, 0, 100, 100); backgroundChanged = true; pixels = null; if (useBufferedImage) { try { /* simulate the following code using reflection: dbimage = new BufferedImage(d.width, d.height, BufferedImage.TYPE_INT_RGB); DataBuffer db = (DataBuffer)(((BufferedImage)dbimage). getRaster().getDataBuffer()); DataBufferInt dbi = (DataBufferInt) db; pixels = dbi.getData(); */ Class biclass = Class.forName("java.awt.image.BufferedImage"); Class dbiclass = Class.forName("java.awt.image.DataBufferInt"); Class rasclass = Class.forName("java.awt.image.Raster"); Constructor cstr = biclass.getConstructor( new Class[] { int.class, int.class, int.class }); backimage = (Image) cstr.newInstance(new Object[] { new Integer(d.width), new Integer(d.height), new Integer(BufferedImage.TYPE_INT_RGB)}); Method m = biclass.getMethod("getRaster", null); Object ras = m.invoke(backimage, null); Object db = rasclass.getMethod("getDataBuffer", null). invoke(ras, null); pixels = (int[]) dbiclass.getMethod("getData", null).invoke(db, null); } catch (Exception ee) { // ee.printStackTrace(); System.out.println("BufferedImage failed"); } } if (pixels == null) { pixels = new int[d.width*d.height]; int i; for (i = 0; i != d.width*d.height; i++) pixels[i] = 0xFF000000; imageSource = new MemoryImageSource(d.width, d.height, pixels, 0, d.width); imageSource.setAnimated(true); imageSource.setFullBufferUpdates(true); backimage = cv.createImage(imageSource); } } void resetDensityGroups() { int i, j, k; for (i = 0; i != densitygridsize; i++) for (j = 0; j != densitygridsize; j++) density[i][j] = 0; int pcount = getParticleCount(); for (i = 0; i != pcount; i++) { Particle p = particles[i]; addToDensityGroup(p); } // invalidate all unused particles so they will be // repositioned if the particle slider is moved to the right, // rather than just falling in where they happened to be last time for (; i != maxParticleCount; i++) { Particle p = particles[i]; p.lifetime = -100; } } int addToDensityGroup(Particle p) { int a = (int)((p.pos[0]+1)*(densitygridsize/2)); int b = (int)((p.pos[1]+1)*(densitygridsize/2)); int n = 0; try { n = ++density[a][b]; if (n > maxParticleCount) System.out.print(a + " " + b + " " + density[a][b] + "\n"); } catch (Exception e) { System.out.print(p.pos[0] + " " + p.pos[1] + "\n"); e.printStackTrace(); } return n; } void removeFromDensityGroup(Particle p) { int a = (int)((p.pos[0]+1)*(densitygridsize/2)); int b = (int)((p.pos[1]+1)*(densitygridsize/2)); try { if (--density[a][b] < 0) System.out.print(a + " " + b + " " + density[a][b] + "\n"); } catch (Exception e) { System.out.print(p.pos[0] + " " + p.pos[1] + "\n"); e.printStackTrace(); } } void positionParticle(Particle p) { int x, y; int bestx = 0, besty = 0; int best = 10000; // we avoid scanning the grid in the same order every time // so that we treat equal-density squares as equally as possible. int randaddx = getrand(densitygridsize); int randaddy = getrand(densitygridsize); for (x = 0; x != densitygridsize; x++) for (y = 0; y != densitygridsize; y++) { int ix = (randaddx + x) % densitygridsize; int iy = (randaddy + y) % densitygridsize; if (density[ix][iy] <= best) { bestx = ix; besty = iy; best = density[ix][iy]; } } p.pos[0] = bestx*densitygroupsize + getrand(100)*densitygroupsize/100.0 - 1; p.pos[1] = besty*densitygroupsize + getrand(100)*densitygroupsize/100.0 - 1; p.lifetime = curfunc.redistribute() ? 500 : 5000; p.stepsize = 1; p.theta = (getrand(101)-50)*pi/50.; p.phi = (getrand(101)-50)*pi/50.; int j; for (j = 0; j != 3; j++) p.vel[j] = 0; } int getParticleCount() { return partCountBar.getValue(); } void resetParticles() { int pcount = getParticleCount(); int i, j; for (i = 0; i != pcount; i++) { Particle p = particles[i]; for (j = 0; j != 2; j++) { p.pos[j] = getrand(200)/100.0 - 1; p.vel[j] = 0; } p.pos[2] = 0; p.lifetime = i*2; p.stepsize = 1; } integralX = -1; resetDensityGroups(); } void kickParticles() { int i, j; for (i = 0; i != getParticleCount(); i++) { Particle p = particles[i]; for (j = 0; j != 2; j++) p.vel[j] += (getrand(100)/99.0 - .5) * .04; } } void generateFunction() { int x, y; if (grid == null) grid = new GridElement[gridsize+1][gridsize+1]; curfunc.setupFrame(); divOffset = curfunc.getDivOffset(); divRange = curfunc.getDivRange(); double mu, xx, xx2, yy, yy2, r, r1, r2, r3, r4; double levelheight = curfunc.getLevelHeight(); //System.out.print("gf1\n"); for (x = 0; x != gridsize+1; x++) for (y = 0; y != gridsize+1; y++) { GridElement ge = grid[x][y] = new GridElement(); ge.curl = ge.div = ge.height = 0; curfunc.setGrid(ge, x, y); /*ge.normdot = (ge.vecX+ge.vecY+1)*(1/1.73)/ java.lang.Math.sqrt(ge.vecX*ge.vecX+ge.vecY*ge.vecY+1);*/ } curfunc.calcDivergence(); double zval = 2.0/gridsize; for (y = 0; y != gridsize; y++) for (x = 0; x != gridsize; x++) { GridElement ge = grid[x][y]; double vecx = grid[x+1][y].height-ge.height; double vecy = grid[x][y+1].height-ge.height; // light vector = (1,1,1)/1.73 // norm vector = (-dzdx, -dzdy, 1) ge.normdot = (vecx+vecy+zval)*(1/1.73)/ java.lang.Math.sqrt(vecx*vecx+vecy*vecy+zval*zval); } for (x = 0; x != gridsize+1; x++) { grid[gridsize][x] = grid[gridsize-1][x]; grid[x][gridsize] = grid[x][gridsize-1]; } //System.out.print("gf2\n"); functionChanged = false; backgroundChanged = true; } double divOffset, divRange; int computeColor(GridElement ge, double c) { if (c < 0) c = 0; if (c > 1) c = 1; //System.out.print(c + "\n"); c = .5 + c * .5; double value = 0; double range = 10; double offset = 4; switch (floorColorChooser.getSelectedIndex()) { case FC_FIELD: value = ge.vecX*ge.vecX+ge.vecY*ge.vecY; offset = 10; range = 16; if (!ge.valid) return 0xFF000080; break; case FC_POTENTIAL: value = ge.height-curfunc.getLevelHeight(); offset = 1; range = 2; break; case FC_CURL: value = ge.curl; offset = 4; range = 10; break; case FC_DIV: value = ge.div; offset = divOffset; range = divRange; break; case FC_NONE: if (!ge.valid) return 0xFF000080; break; } value *= floorBrightMult; double redness = (value < 0) ? (java.lang.Math.log(-value)+offset)/range : 0; double grnness = (value > 0) ? (java.lang.Math.log( value)+offset)/range : 0; if (redness > 1) redness = 1; if (grnness > 1) grnness = 1; if (grnness < 0) grnness = 0; if (redness < 0) redness = 0; /*if (redness > 0) System.out.print(redness + "\n");*/ double grayness = (1-(redness+grnness))*c; double gray = .6; int r = (int) ((c*redness+gray*grayness)*255); int g = (int) ((c*grnness+gray*grayness)*255); int b = (int) ((gray*grayness)*255); return (255<<24) | (r<<16) | (g<<8) | b; } void reinit() { handleResize(); resetParticles(); functionChanged = backgroundChanged = true; } void centerString(Graphics g, String s, int y) { FontMetrics fm = g.getFontMetrics(); g.drawString(s, (winSize.width-fm.stringWidth(s))/2, y); } int shadowBufferTop[], shadowBufferBottom[], shadowBufferTop2[], shadowBufferBottom2[]; final double floorBrightMult = 2; void drawBackground() { if (isFlat) { int x, y; for (y = 0; y < gridsize; y++) for (x = 0; x < gridsize; x++) { GridElement ge = grid[x][y]; int nx = x*winSize.width/gridsize; int ny = winSize.height-(y+1)*winSize.height/gridsize; int nx1 = (x+1)*winSize.width/gridsize; int ny1 = winSize.height-y*winSize.height/gridsize; int col = computeColor(ge, 0); fillRectangle(nx, ny, nx1, ny1, col); ge.visible = true; } drawFloor(); functionChanged = backgroundChanged = false; if (imageSource != null) imageSource.newPixels(); return; } scaleworld(); int x, y; int xdir, xstart, xend; int ydir, ystart, yend; int sc = gridsize; // figure out what order to render the grid elements so that // the ones in back are rendered first. if (viewAngleCos < 0) { ystart = sc; yend = 0; ydir = -1; } else { ystart = 0; yend = sc; ydir = 1; } if (viewAngleSin < 0) { xstart = 0; xend = sc; xdir = 1; } else { xstart = sc; xend = 0; xdir = -1; } boolean xFirst = (-viewAngleSin * xdir > viewAngleCos * ydir); shadowBufferBottom = new int[winSize.width]; shadowBufferTop = new int[winSize.width]; shadowBufferBottom2 = new int[winSize.width]; shadowBufferTop2 = new int[winSize.width]; for (x = 0; x != winSize.width; x++) { shadowBufferBottom[x] = shadowBufferBottom2[x] = 0; shadowBufferTop[x] = shadowBufferTop2[x] = winSize.height-1; } for (x = 0; x != winSize.width*winSize.height; x++) pixels[x] = 0xFF000000; int goffx = (xdir == 1) ? 0 : -1; int goffy = (ydir == 1) ? 0 : -1; //System.out.print("drawing 1\n"); for (x = xstart; x != xend; x += xdir) { for (y = ystart; y != yend; y += ydir) { if (!xFirst) x = xstart; for (; x != xend; x += xdir) { double nx = x*(2.0/gridsize)-1; double ny = y*(2.0/gridsize)-1; double nx1 = (x+xdir)*(2.0/gridsize)-1; double ny1 = (y+ydir)*(2.0/gridsize)-1; map3d(nx, ny, grid[x] [y].height, xpoints, ypoints, 0); map3d(nx1, ny, grid[x+xdir][y].height, xpoints, ypoints, 1); map3d(nx, ny1, grid[x][y+ydir].height, xpoints, ypoints, 2); map3d(nx1, ny1, grid[x+xdir][y+ydir].height, xpoints, ypoints, 3); GridElement ge = grid[x+goffx][y+goffy]; int col = computeColor(ge, ge.normdot); // if (((x&1)^(y&1)) == 0) col ^= 0xff; fillTriangle(xpoints[0], ypoints[0], xpoints[1], ypoints[1], xpoints[3], ypoints[3], col); fillTriangle(xpoints[0], ypoints[0], xpoints[2], ypoints[2], xpoints[3], ypoints[3], col); int cx = (xpoints[0]+xpoints[3])/2; int cy = (ypoints[0]+ypoints[3])/2; boolean vis = false; if (cx >= 0 && cx < winSize.width && cy <= shadowBufferTop[cx] && cy >= 0) vis = true; ge.visible = vis; if (xFirst) break; } if (!xFirst) { int i; for (i = 0; i != winSize.width; i++) { shadowBufferTop[i] = shadowBufferTop2[i]; shadowBufferBottom[i] = shadowBufferBottom2[i]; } } } if (!xFirst) break; int i; for (i = 0; i != winSize.width; i++) { shadowBufferTop[i] = shadowBufferTop2[i]; shadowBufferBottom[i] = shadowBufferBottom2[i]; } } //System.out.print("drawing 2\n"); drawFloor(); //System.out.print("drawing 3\n"); functionChanged = backgroundChanged = false; if (imageSource != null) imageSource.newPixels(); } void drawFloor() { int x, y; switch (floorLineChooser.getSelectedIndex()) { case FL_NONE: break; case FL_GRID: for (x = 0; x != gridsize; x++) for (y = 0; y != gridsize; y += 10) { double nx = x*(2.0/gridsize)-1; double nx1 = (x+1)*(2.0/gridsize)-1; double ny = y*(2.0/gridsize)-1; if (grid[x][y].visible) { map3d(nx, ny, grid[x] [y].height, xpoints, ypoints, 0); map3d(nx1, ny, grid[x+1][y].height, xpoints, ypoints, 1); drawLine(xpoints[0], ypoints[0], xpoints[1], ypoints[1]); } if (grid[y][x].visible) { map3d(ny, nx, grid[y][x].height, xpoints, ypoints, 0); map3d(ny, nx1, grid[y][x+1].height, xpoints, ypoints, 1); drawLine(xpoints[0], ypoints[0], xpoints[1], ypoints[1]); } } break; case FL_EQUIP: if (!curfunc.nonGradient()) renderEquips(); break; case FL_LINES: genLines(); break; } } void fillTriangle(int x1, int y1, int x2, int y2, int x3, int y3, int col) { if (x1 > x2) { if (x2 > x3) { // x1 > x2 > x3 int ay = interp(x1, y1, x3, y3, x2); fillTriangle1(x3, y3, x2, y2, ay, col); fillTriangle1(x1, y1, x2, y2, ay, col); } else if (x1 > x3) { // x1 > x3 > x2 int ay = interp(x1, y1, x2, y2, x3); fillTriangle1(x2, y2, x3, y3, ay, col); fillTriangle1(x1, y1, x3, y3, ay, col); } else { // x3 > x1 > x2 int ay = interp(x3, y3, x2, y2, x1); fillTriangle1(x2, y2, x1, y1, ay, col); fillTriangle1(x3, y3, x1, y1, ay, col); } } else { if (x1 > x3) { // x2 > x1 > x3 int ay = interp(x2, y2, x3, y3, x1); fillTriangle1(x3, y3, x1, y1, ay, col); fillTriangle1(x2, y2, x1, y1, ay, col); } else if (x2 > x3) { // x2 > x3 > x1 int ay = interp(x2, y2, x1, y1, x3); fillTriangle1(x1, y1, x3, y3, ay, col); fillTriangle1(x2, y2, x3, y3, ay, col); } else { // x3 > x2 > x1 int ay = interp(x3, y3, x1, y1, x2); fillTriangle1(x1, y1, x2, y2, ay, col); fillTriangle1(x3, y3, x2, y2, ay, col); } } } int interp(int x1, int y1, int x2, int y2, int x) { if (x1 == x2) return y1; if (x < x1 && x < x2 || x > x1 && x > x2) System.out.print("interp out of bounds\n"); return (int) (y1+((double) x-x1)*(y2-y1)/(x2-x1)); } void fillTriangle1(int x1, int y1, int x2, int y2, int y3, int col) { // x2 == x3 int dir = (x1 > x2) ? -1 : 1; int x = x1; if (x < 0) { x = 0; if (x2 < 0) return; } if (x >= winSize.width) { x = winSize.width-1; if (x2 >= winSize.width) return; } if (y2 > y3) { int q = y2; y2 = y3; y3 = q; } // y2 < y3 while (x != x2+dir) { // XXX this could be speeded up int ya = interp(x1, y1, x2, y2, x); int yb = interp(x1, y1, x2, y3, x); if (ya < 0) ya = 0; if (yb >= winSize.height) yb = winSize.height-1; if (shadowBufferTop2[x] > ya) shadowBufferTop2[x] = ya; if (shadowBufferBottom2[x] < yb) shadowBufferBottom2[x] = yb; int sb1 = shadowBufferTop[x]; int sb2 = shadowBufferBottom[x]; if (!(ya >= sb1 && yb <= sb2)) { for (; ya <= yb; ya++) { if (ya < sb1 || ya > sb2) pixels[x+ya*winSize.width] = col; } } x += dir; if (x < 0 || x >= winSize.width) return; } } void fillRectangle(int x1, int y1, int x2, int y2, int col) { int x, y; for (y = y1; y < y2; y++) for (x = x1; x < x2; x++) pixels[x+y*winSize.width] = col; } void drawLine(int x1, int y1, int x2, int y2) { if (x1 == x2 && y1 == y2) return; if (abs(y2-y1) > abs(x2-x1)) { // y difference is greater, so we step along y's // from min to max y and calculate x for each step int sgn = sign(y2-y1); int x, y; for (y = y1; y != y2+sgn; y += sgn) { x = x1+(x2-x1)*(y-y1)/(y2-y1); if (x >= 0 && y >= 0 && x < winSize.width && y < winSize.height) pixels[x+y*winSize.width] = 0xFFC0C0C0; } } else { // x difference is greater, so we step along x's // from min to max x and calculate y for each step int sgn = sign(x2-x1); int x, y; for (x = x1; x != x2+sgn; x += sgn) { y = y1+(y2-y1)*(x-x1)/(x2-x1); if (x >= 0 && y >= 0 && x < winSize.width && y < winSize.height) pixels[x+y*winSize.width] = 0xFFC0C0C0; } } } int abs(int x) { return x < 0 ? -x : x; } int sign(int x) { return (x < 0) ? -1 : (x == 0) ? 0 : 1; } int min(int a, int b) { return (a < b) ? a : b; } int max(int a, int b) { return (a > b) ? a : b; } double min(double a, double b) { return (a < b) ? a : b; } double max(double a, double b) { return (a > b) ? a : b; } // render equipotentials void renderEquips() { int x, y; for (x = 0; x != gridsize; x++) for (y = 0; y != gridsize; y++) { if (!grid[x][y].visible) continue; // try all possible edge combinations tryEdge(x, y, x+1, y, x, y+1, x+1, y+1); tryEdge(x, y, x+1, y, x, y, x, y+1); tryEdge(x, y, x+1, y, x+1, y, x+1, y+1); tryEdge(x, y, x, y+1, x+1, y, x+1, y+1); tryEdge(x, y, x, y+1, x, y+1, x+1, y+1); tryEdge(x+1, y, x+1, y+1, x, y+1, x+1, y+1); } } class FloatPair { public double x, y; } // interpolate between two points void interpPoint(GridElement ep1, GridElement ep2, int x1, int y1, int x2, int y2, double pval, FloatPair pos) { double interp2 = (pval-ep1.height)/(ep2.height-ep1.height); double interp1 = 1-interp2; pos.x = (x1*interp1 + x2*interp2) * 2./gridsize - 1; pos.y = (y1*interp1 + y2*interp2) * 2./gridsize - 1; } // check to see if pval is between the potential at ep1 and ep2 boolean spanning(GridElement ep1, GridElement ep2, double pval) { if (ep1.height == ep2.height) return false; return !((ep1.height < pval && ep2.height < pval) || (ep1.height > pval && ep2.height > pval)); } // try to draw any equipotentials between edge (x1,y1)-(x2,y2) // and edge (x3,y3)-(x4,y4). void tryEdge(int x1, int y1, int x2, int y2, int x3, int y3, int x4, int y4) { int i; double emult = 5; double mult = 1/(40 * emult * .1); GridElement ep1 = grid[x1][y1]; GridElement ep2 = grid[x2][y2]; GridElement ep3 = grid[x3][y3]; GridElement ep4 = grid[x4][y4]; double pmin = min(min(ep1.height, ep2.height), min(ep3.height, ep4.height)); double pmax = max(max(ep1.height, ep2.height), max(ep3.height, ep4.height)); if (pmin < -5) // was 10 pmin = -5; if (pmax > 5) pmax = 5; int imin = (int) (pmin/mult); int imax = (int) (pmax/mult); for (i = imin; i <= imax; i++) { double pval = i*mult; if (!(spanning(ep1, ep2, pval) && spanning(ep3, ep4, pval))) continue; FloatPair pa = new FloatPair(); FloatPair pb = new FloatPair(); interpPoint(ep1, ep2, x1, y1, x2, y2, pval, pa); interpPoint(ep3, ep4, x3, y3, x4, y4, pval, pb); map3d(pa.x, pa.y, pval, xpoints, ypoints, 0); map3d(pb.x, pb.y, pval, xpoints, ypoints, 1); drawLine(xpoints[0], ypoints[0], xpoints[1], ypoints[1]); } } void drawLineBackground(Graphics g) { int x, y; for (x = 0; x != gridsize; x++) for (y = 0; y != gridsize; y++) grid[x][y].visible = true; if (isFlat) return; for (y = gridsize-1; y >= 0; y--) { for (x = 0; x < gridsize; x += 5) { double ny = y*(2.0/gridsize)-1; double nx1 = (x+1)*(2.0/gridsize)-1; double ny1 = (y+1)*(2.0/gridsize)-1; map3d(nx1, ny, grid[x][y].height, xpoints, ypoints, 1); map3d(nx1, ny1, grid[x][y+1].height, xpoints, ypoints, 2); ypoints[1] = bound_y(ypoints[1]); ypoints[2] = bound_y(ypoints[2]); g.drawLine(xpoints[1], ypoints[1], xpoints[2], ypoints[2]); } } for (y = 0; y < gridsize; y += 5) { for (x = gridsize-1; x >= 0; x--) { double nx = x*(2.0/gridsize)-1; double nx1 = (x+1)*(2.0/gridsize)-1; double ny1 = (y+1)*(2.0/gridsize)-1; map3d(nx, ny1, grid[x] [y].height, xpoints, ypoints, 3); map3d(nx1, ny1, grid[x+1][y].height, xpoints, ypoints, 2); ypoints[3] = bound_y(ypoints[3]); ypoints[2] = bound_y(ypoints[2]); g.drawLine(xpoints[3], ypoints[3], xpoints[2], ypoints[2]); } } } int bound_y(int y) { if (y < -100) y = -100; if (y > winSize.height+100) y = winSize.height+100; return y; } public void paint(Graphics g) { cv.repaint(); } static final double root2 = 1.4142135623730950488016887242096981; double scalex, scaley; void map3d(double x, double y, double z, int xpoints[], int ypoints[], int pt) { map3d(x, y, z, xpoints, ypoints, pt, viewMain); } void map3d(double x, double y, double z, int xpoints[], int ypoints[], int pt, Rectangle view) { if (isFlat) { xpoints[pt] = view.x+(int) ((x+1)*view.width/2); ypoints[pt] = view.y+(int) ((1-y)*view.height/2); return; } if (z < -1000) z = -1000; if (z > 1000) z = 1000; /*if (!(z > -2000 && z < 2000)) System.out.print(z + " in map3d\n");*/ double realx = x*viewAngleCos + y*viewAngleSin; // range: [-10,10] double realy = z-viewHeight; double realz = y*viewAngleCos - x*viewAngleSin + viewDistance; scalex = viewZoom * (view.width/4) * viewDistance; scaley = scalex; int yoff = (int) (scaley*(viewHeight-curfunc.getLevelHeight())/viewDistance); xpoints[pt] = view.x+view.width/2 + (int) (scalex*realx/realz); ypoints[pt] = view.y+view.height/2 - yoff - (int) (scaley*realy/realz); } void scaleworld() { // scalex = viewZoom * (winSize.width/4) * viewDistance; // scaley = scalex; } double getHeight(double x, double y) { x = (x+1)*(gridsize/2); y = (y+1)*(gridsize/2); int ix = (int) x; int iy = (int) y; if (ix >= gridsize || iy >= gridsize) return grid[ix][iy].height; double fracx = x-ix; double fracy = y-iy; return grid[ix][iy].height * (1-fracx) * (1-fracy) + grid[ix+1][iy].height * fracx * (1-fracy) + grid[ix][iy+1].height * (1-fracx) * fracy + grid[ix+1][iy+1].height * fracx * fracy; } void sayCalculating(Graphics realg) { realg.setColor(cv.getBackground()); FontMetrics fm = realg.getFontMetrics(); String s = "Calculating..."; realg.fillRect(0, winSize.height-30, 20+fm.stringWidth(s), 30); realg.setColor(Color.white); realg.drawString(s, 10, winSize.height-10); } long lastTime; double timeStep, partMult; boolean slowDragView = true; public void updateVecDemo(Graphics realg) { Graphics g = dbimage.getGraphics(); if (winSize == null || winSize.width == 0) return; if (xpoints == null) return; checkFlatState(); barFieldStrength = fieldStrength = java.lang.Math.exp((strengthBar.getValue()-50)/10.); if (functionChanged || backgroundChanged) { if (functionChanged) { // && !(draggingView && isFlat)) { sayCalculating(realg); generateFunction(); } if (!slowDragView || !draggingView) { long tm1 = System.currentTimeMillis(); sayCalculating(realg); drawBackground(); long tm2 = System.currentTimeMillis(); slowDragView = (tm2-tm1 > 40); } } scaleworld(); if ((draggingView && slowDragView) || functionChanged) { g.setColor(isFlat ? fieldColors[0] : cv.getBackground()); g.fillRect(0, 0, winSize.width, winSize.height); g.setColor(cv.getForeground()); drawLineBackground(g); } else g.drawImage(backimage, 0, 0, this); boolean allquiet = true; curfunc.setupFrame(); fieldStrength = barFieldStrength; partMult = fieldStrength * reverse * timeStep; int disp = dispChooser.getSelectedIndex(); timeStep = 1; if (!stoppedCheck.getState()) { if (lastTime > 0) timeStep = (System.currentTimeMillis()-lastTime)*.03; if (timeStep > 3) timeStep = 3; lastTime = System.currentTimeMillis(); if (disp != DISP_VECTORS && disp != DISP_NONE) { moveParticles(); allquiet = false; } currentStep += reverse; if (currentStep < 0) currentStep += 800; } else lastTime = 0; if (disp == DISP_VECTORS) drawVectors(g); else if (disp != DISP_NONE) drawParticles(g); g.setColor(Color.gray); if (!isFlat) drawAxes(g); curfunc.finishFrame(); int mode = modeChooser.getSelectedIndex(); if (mode == MODE_LINE_INT) lineIntegral(g, true); else if (mode == MODE_SURF_INT) lineIntegral(g, false); if (parseError) centerString(g, "Can't parse expression", winSize.height-20); //g.drawString("FR: " + framerate, 0, 10); realg.drawImage(dbimage, 0, 0, this); long t = System.currentTimeMillis(); frames++; if (firsttime == 0) firsttime = t; else if (t-firsttime > 1000) { framerate = frames; firsttime = t; frames = 0; } if (!stoppedCheck.getState() && !allquiet) cv.repaint(pause); } void drawAxes(Graphics g) { g.setColor(Color.white); map3d(0, 0, 0, xpoints, ypoints, 0, viewAxes); map3d(1, 0, 0, xpoints, ypoints, 1, viewAxes); drawArrow(g, "x", xpoints[0], ypoints[0], xpoints[1], ypoints[1]); map3d(0, 1, 0, xpoints, ypoints, 1, viewAxes); drawArrow(g, "y", xpoints[0], ypoints[0], xpoints[1], ypoints[1]); map3d(0, 0, 1, xpoints, ypoints, 1, viewAxes); drawArrow(g, "z", xpoints[0], ypoints[0], xpoints[1], ypoints[1]); } void drawVectors(Graphics g) { int x, y, z; DrawData dd = new DrawData(); dd.mult = barFieldStrength*40; dd.g = g; dd.field = new double[3]; dd.vv = new double[3]; vectorSpacing = vecDensityBar.getValue(); double vec[] = new double[3]; vecCount = 0; for (x = 0; x != vectorSpacing; x++) { vec[0] = x*(2.0/(vectorSpacing-1))-1; for (y = 0; y != vectorSpacing; y++) { vec[1] = y*(2.0/(vectorSpacing-1))-1; drawVector(dd, vec); } } } void lineIntegral(Graphics g, boolean line) { if (integralX == -1) return; if (dragStartX == integralX || dragStartY == integralY) return; int x1 = min(dragStartX, integralX); int y1 = min(dragStartY, integralY); int x2 = max(dragStartX, integralX); int y2 = max(dragStartY, integralY); int step = 15; int x; double pos[] = rk_k2; if (!line) { g.setColor(Color.white); g.drawRect(x1, y1, x2-x1+1, y2-y1+1); } double y1p = 1 - 2.*y1/winSize.height; double y2p = 1 - 2.*y2/winSize.height; for (x = x1; x <= x2; x += step) { int step1 = x2-x; if (step1 > step) step1 = step; pos[0] = 2.*x/winSize.width - 1; pos[1] = y1p; lineIntegralStep(g, x, y1, pos, step1, 0, line); pos[1] = y2p; lineIntegralStep(g, x+step1, y2, pos, -step1, 0, line); } int y; double x1p = 2.*x1/winSize.width - 1; double x2p = 2.*x2/winSize.width - 1; for (y = y2; y >= y1; y -= step) { int step1 = y-y1; if (step1 > step) step1 = step; pos[0] = x1p; pos[1] = 1 - 2.*y/winSize.height; lineIntegralStep(g, x1, y, pos, 0, step1, line); pos[0] = x2p; lineIntegralStep(g, x2, y-step1, pos, 0, -step1, line); } boundCheck = false; pos[1] = y1p; double iv1 = numIntegrate(pos, 0, x1p, x2p, line); pos[1] = y2p; double iv2 = numIntegrate(pos, 0, x1p, x2p, line); pos[0] = x1p; double iv3 = numIntegrate(pos, 1, y1p, y2p, line); pos[0] = x2p; double iv4 = numIntegrate(pos, 1, y1p, y2p, line); double ivtot = -iv1+iv2+iv3-iv4; NumberFormat nf = NumberFormat.getInstance(); nf.setMaximumFractionDigits(3); if (ivtot < 1e-7 && ivtot > -1e-7) ivtot = 0; ivtot *= reverse; String s = ((!line) ? "Flux = " : "Circulation = "); s += nf.format(ivtot*1e5); g.setColor(cv.getBackground()); FontMetrics fm = g.getFontMetrics(); g.fillRect(0, winSize.height-30, 20+fm.stringWidth(s), 30); g.setColor(Color.white); g.drawString(s, 10, winSize.height-10); } double numIntegrate(double pos[], int n1, double x1, double x2, boolean line) { int steps = 8; double lastres = 0; double res = 0; int n2 = (line) ? n1 : 1-n1; while (true) { int i; double h = (x2-x1)/steps; res = 0; for (i = 0; i <= steps; i++) { pos[n1] = x1+i*h; double field[] = rk_k1; curfunc.getField(field, pos); int ss = (i == 0 || i == steps) ? 1 : ((i & 1) == 1) ? 4 : 2; res += field[n2]*h*ss; } res /= 3; if (java.lang.Math.abs(lastres-res) < 1e-7) break; lastres = res; steps *= 2; if (steps == 65536) break; } if (!line && n1 == 0) res = -res; return res; } void lineIntegralStep(Graphics g, int x, int y, double pos[], int dx, int dy, boolean line) { double field[] = rk_k1; curfunc.getField(field, pos); double f = (line) ? field[0]*dx+field[1]*dy : field[0]*dy - field[1]*dx; f *= reverse; double dn = java.lang.Math.abs(f*100); if (dn > 1) dn = 1; int col1 = (int) (dn*128+127); int col2 = (int) (127-dn*127); if (!line) { x += dx/2; y -= dy/2; } if (f == 0) { g.setColor(new Color(col2, col2, col2)); g.drawLine(x, y, x+dx, y-dy); } else if (f > 0) { g.setColor(new Color(col1, col2, col2)); if (line) drawArrow(g, null, x, y, x+dx, y-dy); else drawArrow(g, null, x, y, x+dy, y+dx); } else { g.setColor(new Color(col2, col1, col2)); if (line) drawArrow(g, null, x+dx, y-dy, x, y); else drawArrow(g, null, x, y, x-dy, y-dx); } } void genLines() { int i; int lineGridSize = 8; // lineDensityBar.getValue(); if (lineGridSize < 3) lineGridSize = 3; if (lineGridSize > 8) lineGridSize = 8; lineGridSize *= 2; int ct = 30*lineGridSize*lineGridSize; double brightmult = 80*barFieldStrength; fieldStrength = 10; boolean lineGrid[][] = new boolean[lineGridSize][lineGridSize]; double lineGridMult = lineGridSize/2.; double origp[] = new double[3]; double field[] = new double[3]; Particle p = new Particle(); p.lifetime = -1; p.stepsize = 10; int dir = -1; int segs = 0; double lastdist = 0; for (i = 0; i != ct; i++) { if (p.lifetime < 0) { p.lifetime = 1; p.stepsize = 10; segs = 0; lastdist = 0; if (dir == 1) { int j; for (j = 0; j != 3; j++) p.pos[j] = origp[j]; dir = -1; continue; } dir = 1; int px = 0, py = 0; while (true) { if (!lineGrid[px][py]) break; if (++px < lineGridSize) continue; px = 0; if (++py < lineGridSize) continue; break; } if (py == lineGridSize) break; lineGrid[px][py] = true; double offs = .5/lineGridMult; origp[0] = p.pos[0] = px/lineGridMult-1+offs; origp[1] = p.pos[1] = py/lineGridMult-1+offs; } double p1x = p.pos[0]; double p1y = p.pos[1]; double p1z = getHeight(p1x, p1y); GridElement ge = grid[(int) ((p1x+1)*gridsize/2)] [(int) ((p1y+1)*gridsize/2)]; if (!ge.visible) { p.lifetime = -1; continue; } double x[] = p.pos; lineSegment(p, dir); //System.out.print(x[0] + " " + x[1] + " " + x[2] + "\n"); if (p.lifetime < 0) continue; int gx = (int) ((x[0]+1)*lineGridMult); int gy = (int) ((x[1]+1)*lineGridMult); if (!lineGrid[gx][gy]) segs--; lineGrid[gx][gy] = true; ge = grid[(int) ((p.pos[0]+1)*gridsize/2)] [(int) ((p.pos[1]+1)*gridsize/2)]; if (!ge.visible) { p.lifetime = -1; continue; } double dn = brightmult*p.phi; if (dn > 2) dn = 2; // fv.col = (int) (dn*255); map3d(p1x, p1y, p1z, xpoints, ypoints, 0); map3d(p.pos[0], p.pos[1], getHeight(p.pos[0], p.pos[1]), xpoints, ypoints, 1); drawLine(xpoints[0], ypoints[0], xpoints[1], ypoints[1]); double d2 = dist2(origp, x); if (d2 > lastdist) lastdist = d2; else segs++; if (segs > 10 || d2 < .001) p.lifetime = -1; } //System.out.print("vc " + vecCount + " " + ct + "\n"); } // draw the appropriate field vector at xx,yy,zz void drawVector(DrawData dd, double vec[]) { double field[] = dd.field; // calculate field vector curfunc.getField(field, vec); double dn = java.lang.Math.sqrt(field[0]*field[0]+field[1]*field[1]); double dnr = dn*reverse; if (dn > 0) { field[0] /= dnr; field[1] /= dnr; } dn *= dd.mult; if (dn > 2) dn = 2; int col = (int) (dn*255); double sw2 = 1./(vectorSpacing-1); map3d(vec[0], vec[1], 0, xpoints, ypoints, 0); map3d(vec[0]+sw2*field[0], vec[1]+sw2*field[1], 0, xpoints, ypoints, 1); dd.g.setColor(fieldColors[col]); drawArrow(dd.g, null, xpoints[0], ypoints[0], xpoints[1], ypoints[1], 2); } void moveParticles() { int bestd = 0; int i; int pcount = getParticleCount(); for (i = 0; i != pcount; i++) { Particle pt = particles[i]; removeFromDensityGroup(pt); moveParticle(pt); double x[] = pt.pos; if (!(x[0] >= -1 && x[0] < 1 && x[1] >= -1 && x[1] < 1) || (pt.lifetime -= timeStep) < 0) positionParticle(pt); int d = addToDensityGroup(pt); if (d > bestd) bestd = d; } boolean withforce = (dispChooser.getSelectedIndex() == DISP_PART_FORCE); int maxd = (6*getParticleCount()/(densitygridsize*densitygridsize)); if (!withforce && curfunc.redistribute() && bestd > maxd) redistribute(bestd); } static int frames = 0; static int framerate = 0; static long firsttime = 0; void drawParticles(Graphics g) { g.setColor(Color.white); int disp = dispChooser.getSelectedIndex(); if (disp == DISP_VECTORS) { int i; for (i = 0; i != vecCount; i++) { FieldVector fv = vectors[i]; g.setColor(fieldColors[fv.col]); drawArrow(g, null, fv.sx1, fv.sy1, fv.sx2, fv.sy2, 2); } return; } int pcount = getParticleCount(); int i; wooft += .3; if (disp == DISP_CURLERS) pcount = (pcount+4)/5; for (i = 0; i < pcount; i++) { Particle p = particles[i]; double pos[] = p.pos; GridElement ge = grid[(int) ((pos[0]+1)*gridsize/2)] [(int) ((pos[1]+1)*gridsize/2)]; map3d(pos[0], pos[1], getHeight(pos[0], pos[1]), xpoints, ypoints, 0); if (xpoints[0] < 0 || xpoints[0] >= winSize.width || ypoints[0] < 0 || ypoints[0] >= winSize.height) continue; if (disp == DISP_CURLERS) { g.setColor(p.color); final double len = .02; double ax = java.lang.Math.cos(p.theta)*len; double ay = java.lang.Math.sin(p.theta)*len; double offx = ax; double offy = ay; double a1 = curlcalc(p.pos[0]+offx, p.pos[1]+offy, -ay, ax); double a2 = curlcalc(p.pos[0]-offy, p.pos[1]+offx, -ax, -ay); double a3 = curlcalc(p.pos[0]-offx, p.pos[1]-offy, ay, -ax); double a4 = curlcalc(p.pos[0]+offy, p.pos[1]-offx, ax, ay); p.theta += (a1+a2+a3+a4)/(4*len*len); map3d(p.pos[0]-offx, p.pos[1]-offy, 0, xpoints, ypoints, 0); map3d(p.pos[0]+offx, p.pos[1]+offy, 0, xpoints, ypoints, 1); map3d(p.pos[0]-offy, p.pos[1]+offx, 0, xpoints, ypoints, 2); map3d(p.pos[0]+offy, p.pos[1]-offx, 0, xpoints, ypoints, 3); g.drawLine(xpoints[0], ypoints[0], xpoints[1], ypoints[1]); g.drawLine(xpoints[2], ypoints[2], xpoints[3], ypoints[3]); g.fillOval(xpoints[0]-1, ypoints[0]-1, 3, 3); } else if (ge.visible && ge.valid) g.fillRect(xpoints[0], ypoints[0]-1, 2, 2); } } double wooft = 0; void drawPlane(Graphics g, double sizex, double sizey, double z) { g.setColor(darkYellow); map3d(-sizex, -sizey, z, xpoints, ypoints, 0); map3d(+sizex, -sizey, z, xpoints, ypoints, 1); map3d(+sizex, +sizey, z, xpoints, ypoints, 2); map3d(-sizex, +sizey, z, xpoints, ypoints, 3); g.fillPolygon(xpoints, ypoints, 4); } void drawArrow(Graphics g, String text, int x1, int y1, int x2, int y2) { drawArrow(g, text, x1, y1, x2, y2, 5); } void drawArrow(Graphics g, String text, int x1, int y1, int x2, int y2, int as) { g.drawLine(x1, y1, x2, y2); double l = java.lang.Math.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); if (l > as/2) { // was as double hatx = (x2-x1)/l; double haty = (y2-y1)/l; g.drawLine(x2, y2, (int) (haty*as-hatx*as+x2), (int) (-hatx*as-haty*as+y2)); g.drawLine(x2, y2, (int) (-haty*as-hatx*as+x2), (int) (hatx*as-haty*as+y2)); if (text != null) g.drawString(text, (int) (x2+hatx*10), (int) (y2+haty*10)); } } int rediscount; void redistribute(int mostd) { if (mostd < 5) return; rediscount++; int maxd = (6*getParticleCount()/ (densitygridsize*densitygridsize)); int i; int pn = 0; int pcount = getParticleCount(); for (i = rediscount % 4; i < pcount; i+=4) { Particle p = particles[i]; int a = (int)((p.pos[0]+1)*(densitygridsize/2)); int b = (int)((p.pos[1]+1)*(densitygridsize/2)); if (density[a][b] <= maxd) continue; p.lifetime = -1; pn++; } //System.out.print("redist " + mostd + " " + pn + "\n"); } double curlcalc(double x, double y, double ax, double ay) { rk_yn[0] = x; rk_yn[1] = y; curfunc.getField(rk_k1, rk_yn); return partMult*(rk_k1[0]*ax + rk_k1[1]*ay); } double distance(Particle p) { return distance(p.pos[0], p.pos[1]); } double distance(double y[]) { return distance(y[0], y[1]); } double distance(double x, double y) { return java.lang.Math.sqrt(x*x+y*y+.000000001); } void rotateParticleAdd(double result[], double y[], double mult, double cx, double cy) { result[0] += -mult*(y[1]-cy); result[1] += mult*(y[0]-cx); } void rotateParticle(double result[], double y[], double mult) { result[0] = -mult*y[1]; result[1] = mult*y[0]; result[2] = 0; } void edit(MouseEvent e) { int x = e.getX(); int y = e.getY(); editView(x, y); } void editView(int x, int y) { if (modeChooser.getSelectedIndex() == MODE_VIEW_ROTATE) { if (isFlat) return; viewAngle = (dragStartX-x)/40. + viewAngleDragStart; while (viewAngle < 0) viewAngle += 2*pi; while (viewAngle >= 2*pi) viewAngle -= 2*pi; viewAngleCos = java.lang.Math.cos(viewAngle); viewAngleSin = java.lang.Math.sin(viewAngle); viewHeight = -(dragStartY-y)/10. + viewHeightDragStart; if (viewHeight > 9) viewHeight = 9; if (viewHeight < -9) viewHeight = -9; draggingView = backgroundChanged = true; cv.repaint(pause); return; } if (modeChooser.getSelectedIndex() == MODE_VIEW_ZOOM) { if (isFlat) return; viewZoom = (x-dragStartX)/40. + viewZoomDragStart; if (viewZoom < .1) viewZoom = .1; draggingView = backgroundChanged = true; cv.repaint(pause); return; } if (modeChooser.getSelectedIndex() == MODE_LINE_INT || modeChooser.getSelectedIndex() == MODE_SURF_INT) { integralX = x; integralY = y; cv.repaint(pause); } } public void componentHidden(ComponentEvent e){} public void componentMoved(ComponentEvent e){} public void componentShown(ComponentEvent e) { cv.repaint(pause); } public void componentResized(ComponentEvent e) { handleResize(); cv.repaint(pause); } public void actionPerformed(ActionEvent e) { vectors = null; if (e.getSource() == resetButton) resetParticles(); if (e.getSource() == kickButton) kickParticles(); if (e.getSource() == infoButton) { String s = curfunc.getClass().getName(); try { s = s.substring(s.lastIndexOf('.')+1); applet.getAppletContext().showDocument( new URL(applet.getCodeBase(), "functions.html" + '#' + s), "functionHelp"); } catch (Exception ex) { } } curfunc.actionPerformed(); cv.repaint(pause); } public boolean handleEvent(Event ev) { if (ev.id == Event.WINDOW_DESTROY) { destroyFrame(); return true; } return super.handleEvent(ev); } void destroyFrame() { if (applet == null) dispose(); else applet.destroyFrame(); } public void scrollbarValueChanged(DecentScrollbar ds) { vectors = null; System.out.print(ds.getValue() + "\n"); if (ds == partCountBar) resetDensityGroups(); if (ds == aux1Bar || ds == aux2Bar || ds == aux3Bar) { functionChanged = true; draggingView = true; } cv.repaint(pause); } public void scrollbarFinished(DecentScrollbar ds) { draggingView = false; cv.repaint(pause); } public void mouseDragged(MouseEvent e) { dragging = true; edit(e); } public void mouseMoved(MouseEvent e) { } public void mouseClicked(MouseEvent e) { } public void mouseEntered(MouseEvent e) { } public void mouseExited(MouseEvent e) { } public void mousePressed(MouseEvent e) { if ((e.getModifiers() & MouseEvent.BUTTON1_MASK) == 0) return; dragStartX = e.getX(); dragStartY = e.getY(); viewAngleDragStart = viewAngle; viewHeightDragStart = viewHeight; viewZoomDragStart = viewZoom; dragging = true; edit(e); } public void mouseReleased(MouseEvent e) { if ((e.getModifiers() & MouseEvent.BUTTON1_MASK) == 0) return; dragging = draggingView = false; cv.repaint(pause); } void dispChooserChanged() { int disp = dispChooser.getSelectedIndex(); showA = false; // XXX if (disp == DISP_PART_FORCE) kickButton.enable(); else kickButton.disable(); potentialLabel.hide(); potentialBar.hide(); vecDensityLabel.hide(); vecDensityBar.hide(); /*lineDensityLabel.hide(); lineDensityBar.hide();*/ partCountLabel.hide(); partCountBar.hide(); strengthLabel.show(); strengthBar.show(); if (disp == DISP_VECTORS) { vecDensityLabel.show(); vecDensityBar.show(); } else { partCountLabel.show(); partCountBar.show(); } validate(); resetParticles(); } public void itemStateChanged(ItemEvent e) { vectors = null; cv.repaint(pause); reverse = (reverseCheck.getState()) ? -1 : 1; if (e.getItemSelectable() == dispChooser) { dispChooserChanged(); resetParticles(); } if (e.getItemSelectable() == functionChooser) functionChanged(); if (e.getItemSelectable() == reverseCheck) functionChanged = true; if (e.getItemSelectable() == floorColorChooser || e.getItemSelectable() == floorLineChooser) backgroundChanged = true; } void checkFlatState() { boolean oldFlat = isFlat; int disp = dispChooser.getSelectedIndex(); isFlat = flatCheck.getState() || curfunc.nonGradient() || disp == DISP_VECTORS || disp == DISP_CURLERS; int mode = modeChooser.getSelectedIndex(); if (mode == MODE_LINE_INT || mode == MODE_SURF_INT) isFlat = true; if (isFlat != oldFlat) backgroundChanged = true; } void functionChanged() { reverse = 1; reverseCheck.setState(false); parseError = false; curfunc = (VecFunction) functionList.elementAt(functionChooser.getSelectedIndex()); int i; for (i = 0; i != 3; i++) { auxBars[i].label.hide(); auxBars[i].bar.hide(); } for (i = 0; i != 2; i++) textFields[i].hide(); if (textFieldLabel != null) textFieldLabel.hide(); strengthBar.setValue(80); //setupDispChooser(!curfunc.nonGradient()); curfunc.setup(); validate(); resetParticles(); dispChooserChanged(); functionChanged = true; integralX = -1; } void setupDispChooser(boolean potential) { dispChooser.removeAll(); dispChooser.add("Display: Particles (Vel.)"); dispChooser.add("Display: Particles (Force)"); dispChooser.add("Display: Field Vectors"); dispChooser.add("Display: None"); #ifdef BUILD_V dispChooser.add("Display: Curl Detectors"); #endif } void setupBar(int n, String text, int val) { auxBars[n].label.setText(text); auxBars[n].label.show(); auxBars[n].bar.setValue(val); auxBars[n].bar.show(); } void cross(double res[], double v1[], double v2[]) { res[0] = v1[1]*v2[2] - v1[2]*v2[1]; res[1] = v1[2]*v2[0] - v1[0]*v2[2]; res[2] = v1[0]*v2[1] - v1[1]*v2[0]; } double dot(double v1[], double v2[]) { return v1[0]*v2[0] + v1[1]*v2[1]; } boolean boundCheck; double oldY[]; double rk_k1[] = new double[6]; double rk_k2[] = new double[6]; double rk_k3[] = new double[6]; double rk_k4[] = new double[6]; double rk_yn[] = new double[6]; void rk(int order, double x, double Y[], double stepsize) { int i; // x is not used... if (order == 2) { // velocity-based motion double fmult = stepsize * partMult; for (i = 0; i != order; i++) rk_yn[i] = Y[i]; curfunc.getField(rk_k1, rk_yn); for (i = 0; i != order; i++) rk_yn[i] = (Y[i] + 0.5*fmult*rk_k1[i]); curfunc.getField(rk_k2, rk_yn); for (i = 0; i != order; i++) rk_yn[i] = (Y[i] + 0.5*fmult*rk_k2[i]); curfunc.getField(rk_k3, rk_yn); for (i = 0; i != order; i++) rk_yn[i] = (Y[i] + fmult*rk_k3[i]); curfunc.getField(rk_k4, rk_yn); for (i = 0; i != order; i++) Y[i] = Y[i] + fmult*(rk_k1[i]+2*(rk_k2[i]+rk_k3[i])+rk_k4[i])/6; Y[2] = rk_k4[2]; } else { // force-based (could share more code with above, but this is // called a lot so we want it to be fast) double fmult = stepsize * partMult; for (i = 0; i != order; i++) rk_yn[i] = Y[i]; getForceField(rk_k1, rk_yn, stepsize, fmult); for (i = 0; i != order; i++) rk_yn[i] = (Y[i] + 0.5*rk_k1[i]); getForceField(rk_k2, rk_yn, stepsize, fmult); for (i = 0; i != order; i++) rk_yn[i] = (Y[i] + 0.5*rk_k2[i]); getForceField(rk_k3, rk_yn, stepsize, fmult); for (i = 0; i != order; i++) rk_yn[i] = (Y[i] + rk_k3[i]); getForceField(rk_k4, rk_yn, stepsize, fmult); for (i = 0; i != order; i++) Y[i] = Y[i] + (rk_k1[i]+2*(rk_k2[i]+rk_k3[i])+rk_k4[i])/6; Y[4] = rk_k4[4]; } } void getForceField(double result[], double y[], double stepsize, double fmult) { // get the field vector for the current function curfunc.getField(result, y); // copy height information result[4] = result[2]; // the field data has been written into result[0:1] where it will // directly influence the position (y[0:1]), but we want it to // influence the velocity (y[2:3]), so we move it. In the // process, we multiply by fmult (to get the reverse button // and field strength slider to work). It would be nice if we // could just say getField(result+2, y), but that's java for you. int i; for (i = 0; i != 2; i++) result[i+2] = fmult*result[i]*.1; // here we fill in result[0:1] so that the particle position will // change according to the velocity. for (i = 0; i != 2; i++) result[i] = stepsize*timeStep*rk_yn[i+2]; } double rk_Y[] = new double[6]; double rk_Yhalf[] = new double[6]; double rk_oldY[] = new double[6]; double ls_fieldavg[] = new double[3]; void moveParticle(Particle p) { int disp = dispChooser.getSelectedIndex(); int numIter=0; double maxh=1; double error=0.0, E = .001, localError; boolean useForce = (disp == DISP_PART_FORCE); int order = useForce ? 4 : 2; double Y[] = rk_Y; double Yhalf[] = rk_Yhalf; oldY = rk_oldY; int i; for (i = 0; i != 2; i++) oldY[i] = Y[i] = Yhalf[i] = p.pos[i]; if (useForce) for (i = 0; i != 2; i++) Y[i+2] = Yhalf[i+2] = p.vel[i]; double t = 0; if (!curfunc.useRungeKutta()) { boundCheck = false; curfunc.getField(Yhalf, Y); if (boundCheck && (!useForce || curfunc.checkBoundsWithForce())) { p.pos[0] = -100; return; } double fmult = partMult; if (useForce) { fmult *= .1; for (i = 0; i != 2; i++) { p.vel[i] += fmult*Yhalf[i]; p.pos[i] += p.vel[i]*timeStep; } } else { for (i = 0; i != 2; i++) p.pos[i] += fmult*Yhalf[i]; } p.pos[2] = Yhalf[2]; for (i = 0; i != 2; i++) Y[i] = p.pos[i]; if (curfunc.checkBounds(Y, oldY)) p.pos[0] = -100; return; } boolean adapt = curfunc.useAdaptiveRungeKutta(); double h = (adapt) ? p.stepsize : 1; int steps = 0; double minh = .0001; while (t >= 0 && t < 1) { if (t+h > 1) h = 1-t; boundCheck = false; // estimate one full step rk(order, 0, Y, h); // bail out after one iteration for some slow fields if (!adapt) break; // estimate two half steps rk(order, 0, Yhalf, h*0.5); rk(order, 0, Yhalf, h*0.5); if (boundCheck && (!useForce || curfunc.checkBoundsWithForce())) { p.pos[0] = -100; return; } // estimate the local error localError = java.lang.Math.abs(Y[0] - Yhalf[0]) + java.lang.Math.abs(Y[1] - Yhalf[1]); if (localError > E && h > minh) { //System.out.print(E + " " + t + " " + localError + " " + h + "\n"); h *= 0.75; // decrease the step size if (h < minh) h = minh; for (i = 0; i != order; i++) Y[i] = Yhalf[i] = oldY[i]; continue; } else if (localError < (E * 0.5)) { h *= 1.25; // increase the step size if (h > maxh) h = maxh; } for (i = 0; i != order; i++) oldY[i] = Yhalf[i] = Y[i]; t += h; steps++; } if (boundCheck && (!useForce || curfunc.checkBoundsWithForce())) { p.pos[0] = -100; return; } p.stepsize = h; for (i = 0; i != 3; i++) p.pos[i] = Y[i]; if (useForce) { for (i = 0; i != 2; i++) p.vel[i] = Y[i+2]; p.pos[2] = Y[4]; } } double dist2(double a[], double b[]) { double c0 = a[0]-b[0]; double c1 = a[1]-b[1]; return c0*c0+c1*c1; } void lineSegment(Particle p, int dir) { int numIter=0; double maxh=20; double error=0.0, E = .001, localError; int order = 2; double Y[] = rk_Y; double Yhalf[] = rk_Yhalf; oldY = rk_oldY; int i; for (i = 0; i != 2; i++) oldY[i] = Y[i] = Yhalf[i] = p.pos[i]; double h = p.stepsize; ls_fieldavg[0] = ls_fieldavg[1] = ls_fieldavg[2] = 0; int steps = 0; double minh = .1; double segSize2max = 1./gridsize; double segSize2min = segSize2max/4; double lastd = 0; int avgct = 0; while (true) { boundCheck = false; steps++; if (steps > 100) { System.out.print("maxsteps\n"); p.lifetime = -1; return; } //System.out.print(h + " " + boundCheck + "/\n"); // estimate one full step rk(order, 0, Y, dir*h); // estimate two half steps rk(order, 0, Yhalf, dir*h*0.5); rk(order, 0, Yhalf, dir*h*0.5); //System.out.print(h + " " + boundCheck + "\n"); if (boundCheck) { for (i = 0; i != order; i++) Y[i] = Yhalf[i] = oldY[i]; h /= 2; if (h < minh) { p.lifetime = -1; return; } continue; } if (Y[0] < -1 || Y[0] >= .999 || Y[1] < -1 || Y[1] >= .999) { for (i = 0; i != order; i++) Y[i] = Yhalf[i] = oldY[i]; h /= 2; if (h < minh) { p.lifetime = -1; return; } continue; } // estimate the local error localError = java.lang.Math.abs(Y[0] - Yhalf[0]) + java.lang.Math.abs(Y[1] - Yhalf[1]); if (localError > E && h > minh) { h *= 0.75; // decrease the step size if (h < minh) h = minh; for (i = 0; i != order; i++) Y[i] = Yhalf[i] = oldY[i]; continue; } else if (localError < (E * 0.5)) { h *= 1.25; // increase the step size if (h > maxh) h = maxh; } double d = dist2(p.pos, Y); if (!(d-lastd > 1e-10)) { // we're not getting anywhere! //System.out.print("nga " + d + " " + lastd + "\n"); p.lifetime = -1; return; } if (d > segSize2max) { h /= 2; if (h < minh) { p.lifetime = -1; return; } for (i = 0; i != order; i++) Y[i] = Yhalf[i] = oldY[i]; continue; } ls_fieldavg[0] += rk_k1[0]; ls_fieldavg[1] += rk_k1[1]; avgct++; if (d > segSize2min) break; lastd = d; for (i = 0; i != order; i++) oldY[i] = Yhalf[i] = Y[i]; } //System.out.print(steps + " ss\n"); p.stepsize = h; for (i = 0; i != 3; i++) p.pos[i] = Y[i]; p.phi = java.lang.Math.sqrt(ls_fieldavg[0]*ls_fieldavg[0]+ ls_fieldavg[1]*ls_fieldavg[1])/avgct; } int doubleToGrid(double x) { return (int) ((x+1)*gridsize/2); } double gridToDouble(int x) { return (x*2./gridsize)-1; } abstract class VecFunction { abstract String getName(); abstract VecFunction createNext(); boolean nonGradient() { return false; } boolean useRungeKutta() { return true; } boolean useAdaptiveRungeKutta() { return true; } boolean checkBoundsWithForce() { return true; } boolean checkBounds(double y[], double oldY[]) { return false; } abstract void getField(double result[], double y[]); boolean redistribute() { return true; } void setup() {} void setupFrame() {} void finishFrame() {} void actionPerformed() {} void calcDivergence() {} double getLevelHeight() { return 0; } void setGrid(GridElement ge, int x, int y) { double xx[] = rk_k1; double res[] = rk_k2; double res1[] = rk_k3; xx[0] = (x*2./gridsize)-1; xx[1] = (y*2./gridsize)-1; xx[2] = 0; boundCheck = false; getField(res, xx); ge.vecX = reverse*res[0]*70; ge.vecY = reverse*res[1]*70; ge.height = reverse*res[2]*.625; ge.valid = !boundCheck; /*if (!(ge.height >= -9999 && ge.height <= 9999)) System.out.print(xx[0] + " " + xx[1] + " " + ge.height + "\n"); */ double xorig0 = xx[0]; xx[0] += 1e-8; getField(res1, xx); ge.div = res1[0] - res[0]; ge.curl = res1[1] - res[1]; xx[0] = xorig0; xx[1] += 1e-8; getField(res1, xx); ge.div = (ge.div + res1[1] - res[1])*1e10*reverse; ge.curl = (ge.curl - (res1[0] - res[0]))*1e10*reverse; } double getDivOffset() { return 4; } double getDivRange() { return 11; } }; class InverseRadial extends VecFunction { double lineLen; String getName() { return BUILD_CASE_EMV("charged line", null, "1/r single line"); } void getField(double result[], double y[]) { double r = distance(y[0], y[1]); if (r < lineWidth) boundCheck = true; double r2 = r*r; result[0] = -.0002*y[0]/r2; result[1] = -.0002*y[1]/r2; result[2] = .4*java.lang.Math.log(r+1e-300); } void setup() { lineLen = 1; } VecFunction createNext() { return new InverseRadialDouble(); } }; class InverseRadialDouble extends VecFunction { InverseRadialDouble() { sign = 1; } String getName() { return BUILD_CASE_EMV("line charge double", null, "1/r double lines"); } double sign; void getField(double result[], double y[]) { // we want to make sure sep is a multiple of the grid spacing // so that the divergence gets calculated correctly. double sep = gridToDouble(gridsize/2+ aux1Bar.getValue()*gridsize/200); double xx1 = y[0] - sep; double xx2 = y[0] + sep; double r1 = distance(xx1, y[1]); double r2 = distance(xx2, y[1]); if (r1 < lineWidth || r2 < lineWidth) boundCheck = true; double q = .0002; double r1s = 1/(r1*r1); double r2s = 1/(r2*r2*sign); result[0] = q*(-xx1 *r1s-xx2 *r2s); result[1] = q*(-y[1]*r1s-y[1]*r2s); result[2] = .2*(java.lang.Math.log(r1+1e-20)+ sign*java.lang.Math.log(r2+1e-20)); } void setup() { setupBar(0, "Line Separation", 30); } VecFunction createNext() { return BUILD_CASE_EMV(new InverseRadialDipole(), null, new InverseSquaredRadial()); } }; #ifdef BUILD_E class InverseRadialDipole extends InverseRadialDouble { InverseRadialDipole() { sign = -1; } String getName() { return "dipole lines"; } VecFunction createNext() { return new InverseRadialQuad(); } }; class InverseRadialQuad extends VecFunction { String getName() { return "quad lines"; } void getField(double result[], double y[]) { // we want to make sure sep is a multiple of the grid spacing // so that the divergence gets calculated correctly. double sep = gridToDouble(gridsize/2+ aux1Bar.getValue()*gridsize/200); double xx1 = y[0] + sep; double xx2 = y[0] - sep; double yy1 = y[1] + sep; double yy2 = y[1] - sep; double r1 = distance(xx1, yy1); double r2 = distance(xx2, yy1); double r3 = distance(xx1, yy2); double r4 = distance(xx2, yy2); if (r1 < lineWidth || r2 < lineWidth || r3 < lineWidth || r4 < lineWidth) boundCheck = true; double q = .0003; result[0] = q*(-xx1/(r1*r1)-xx2/(r4*r4) +xx2/(r2*r2)+xx1/(r3*r3)); result[1] = q*(-yy1/(r1*r1)-yy2/(r4*r4) +yy1/(r2*r2)+yy2/(r3*r3)); result[2] = .2*(+java.lang.Math.log(r1+1e-20) -java.lang.Math.log(r2+1e-20) -java.lang.Math.log(r3+1e-20) +java.lang.Math.log(r4+1e-20)); } void setup() { setupBar(0, "Line Separation", 30); } VecFunction createNext() { return new InverseSquaredRadial(); } }; #endif /* BUILD_E */ class InverseSquaredRadial extends VecFunction { String getName() { return BUILD_CASE_EMV("point charge", null, "1/r^2 single"); } void getField(double result[], double y[]) { double r = distance(y); if (r < chargeSize) boundCheck = true; double r3 = r*r*r; double q = .0003/r3; result[0] = -y[0]*q; result[1] = -y[1]*q; result[2] = -.3/r; } static final double chargeSize = .001; VecFunction createNext() { return new InverseSquaredRadialDouble(); } }; class InverseSquaredRadialDouble extends InverseSquaredRadial { String getName() { return BUILD_CASE_EMV("point charge double", null, "1/r^2 double"); } double sign2; void getField(double result[], double y[]) { // we want to make sure sep is a multiple of the grid spacing // so that the divergence gets calculated correctly. double sep = gridToDouble(gridsize/2+ aux1Bar.getValue()*gridsize/200); double xx1 = y[0]-sep; double xx2 = y[0]+sep; double r1 = distance(xx1, y[1]); if (r1 < chargeSize) boundCheck = true; double r2 = distance(xx2, y[1]); if (r2 < chargeSize) boundCheck = true; double q = .0003; double rq1 = q/(r1*r1*r1); double rq2 = q/(r2*r2*r2) * sign2; result[0] = -xx1 *rq1-xx2 *rq2; result[1] = -y[1]*rq1-y[1]*rq2; result[2] = -.05/r1 - .05*sign2/r2; if (sign2 == -1) result[2] *= 2; } void setup() { sign2 = 1; setupBar(0, "Charge Separation", 30); } VecFunction createNext() { return BUILD_CASE_EMV(new InverseSquaredRadialDipole(), null, new InverseRotational()); } }; #ifdef BUILD_E class InverseSquaredRadialDipole extends InverseSquaredRadialDouble { String getName() { return "dipole"; } void setup() { super.setup(); sign2 = -1; } VecFunction createNext() { return new InverseSquaredRadialQuad(); } }; class InverseSquaredRadialQuad extends InverseSquaredRadial { String getName() { return "quadrupole"; } void getField(double result[], double y[]) { // we want to make sure sep is a multiple of the grid spacing // so that the divergence gets calculated correctly. double sep = gridToDouble(gridsize/2+ aux1Bar.getValue()*gridsize/200); double xx1 = y[0] - sep; double xx2 = y[0] + sep; double yy1 = y[1] - sep; double yy2 = y[1] + sep; double r1 = distance(xx1, yy1); double r2 = distance(xx2, yy1); double r3 = distance(xx1, yy2); double r4 = distance(xx2, yy2); if (r1 < chargeSize || r2 < chargeSize || r3 < chargeSize || r4 < chargeSize) boundCheck = true; double q = .0003; double rq1 = q/(r1*r1*r1); double rq2 = q/(r2*r2*r2); double rq3 = q/(r3*r3*r3); double rq4 = q/(r4*r4*r4); result[0] = -xx1*rq1-xx2*rq4+xx2*rq2+xx1*rq3; result[1] = -yy1*rq1-yy2*rq4+yy1*rq2+yy2*rq3; result[2] = .05*(-1/r1+1/r2+1/r3-1/r4); } void setup() { super.setup(); setupBar(0, "Charge Separation", 30); } VecFunction createNext() { return new ConductingPlate(); } }; class ConductingPlate extends VecFunction { String getName() { return "conducting plate"; } Complex z, z2; boolean plate; double a, base; ConductingPlate() { z = new Complex(); z2 = new Complex(); plate = true; } void setupFrame() { a = (aux1Bar.getValue()+1)/100.; z.set(0, 1/a); z.arcsin(); base = z.b; } void getField(double result[], double y[]) { // smythe p89 if (y[1] >= -.02 && y[1] <= .02) { if ((plate && y[0] >= -a && y[0] <= a) || (!plate && (y[0] >= a || y[0] <= -a))) boundCheck = true; } z.set(y[0]/a, y[1]/a); if (y[1] < 0 && plate) z.b = -z.b; z2.set(z); z2.arcsin(); result[2] = (plate) ? z2.b/base-1 : -z2.a*.6; // here we calculate ((1-(z/a)^2)^(-1/2))/a, which is // d/dz arcsin(z/a) z.square(); z.mult(-1); z.add(1); z.pow(-.5); z.mult(1/a); if (plate) { // field = (Im dw/dz, Re dw/dz) result[1] = z.a * -.0007; result[0] = z.b * -.0007; if (y[1] <= 0) result[1] = -result[1]; } else { // field = (Re dw/dz, -Im dw/dz) result[0] = z.a * .0007; result[1] = -z.b * .0007; if (y[1] == 0) result[1] = -result[1]; } } void setup() { setupBar(0, "Plate Size", 60); } double getDivOffset() { return -17.3; } double getDivRange() { return 2.5; } VecFunction createNext() { return new ChargedPlate(); } }; class ChargedPlate extends ConductingPlate { Complex cz; ChargedPlate() { cz = new Complex(); } String getName() { return "charged plate"; } double getDivOffset() { return 4; } double getDivRange() { return 11; } double getPot(double a1, double a2, double y) { cz.set (y, -a1); cz.mult(y, a2); cz.log(); double b1 = cz.b; cz.set (y, a1); cz.mult(y, -a2); cz.log(); double y2 = y*y; if (y2 == 0) y2 = 1e-8; return .3*(2*(a1-a2) + (b1-cz.b)*y + a2*java.lang.Math.log(a2*a2+y2) - a1*java.lang.Math.log(a1*a1+y2)); } void calcDivergence() { double sep = aux2Bar.getValue()/100.; int x, y; for (x = 0; x != gridsize; x++) { double xx = gridToDouble(x); if (xx < -a || xx > a) continue; grid[x][gridsize/2].div = -reverse; } } void getField(double result[], double y[]) { if (y[1] >= -.01 && y[1] <= .01 && (y[0] >= -a && y[0] <= a)) boundCheck = true; double a1 = -a-y[0]; double a2 = a-y[0]; double y2 = y[1]*y[1]; if (y2 == 0) y2 = 1e-8; double q = .0003/a; result[0] = .5*q* java.lang.Math.log((y2+a2*a2)/(y2+a1*a1)); result[1] = q* (java.lang.Math.atan(a1/y[1])- java.lang.Math.atan(a2/y[1])); result[2] = .4*getPot(a1, a2, y[1])/a; } VecFunction createNext() { return new ChargedPlatePair(); } }; class ChargedPlatePair extends ChargedPlate { String getName() { return "charged plate pair"; } boolean useRungeKutta() { return false; } double dipole; ChargedPlatePair() { dipole = 1; } void getField(double result[], double y[]) { double sep = aux2Bar.getValue()/100.; if ((y[1] >= -.01+sep && y[1] <= .01+sep || y[1] >= -.01-sep && y[1] <= .01-sep) && y[0] >= -a && y[0] <= a) boundCheck = true; double a1 = -a-y[0]; double a2 = a-y[0]; double y1 = y[1]-sep; double y12 = y1*y1; if (y12 == 0) y12 = 1e-8; double y2 = y[1]+sep; double y22 = y2*y2; if (y22 == 0) y22 = 1e-8; double q = .0003/a; result[0] = .5*q* (java.lang.Math.log((y12+a2*a2)/(y12+a1*a1)) + dipole* java.lang.Math.log((y22+a2*a2)/(y22+a1*a1))); result[1] = q* (java.lang.Math.atan(a1/y1) -java.lang.Math.atan(a2/y1) +dipole*(java.lang.Math.atan(a1/y2) -java.lang.Math.atan(a2/y2))); result[2] = .4*(getPot(a1, a2, y1)+dipole*getPot(a1, a2, y2))/a; } void calcDivergence() { double sep = aux2Bar.getValue()/100.; int x, y; for (x = 0; x != gridsize; x++) { double xx = gridToDouble(x); if (xx < -a || xx > a) continue; y = doubleToGrid(sep); grid[x][y].div = -reverse; y = doubleToGrid(-sep); grid[x][y].div = -dipole*reverse; } } void setup() { setupBar(0, "Sheet Size", 60); setupBar(1, "Sheet Separation", 33); } boolean checkBounds(double y[], double oldY[]) { double size = aux1Bar.getValue()/100.; double sep = aux2Bar.getValue()/100.; // check to see if particle has passed through one of plates if (y[0] >= -size && y[0] <= size) { if (y[1] > sep) { if (oldY[1] < sep) return true; } else if (y[1] < -sep) { if (oldY[1] > -sep) return true; } else if (oldY[1] > sep || oldY[1] < -sep) return true; } return false; } VecFunction createNext() { return new ChargedPlateDipole(); } }; class ChargedPlateDipole extends ChargedPlatePair { String getName() { return "charged plate dipole"; } ChargedPlateDipole() { dipole = -1; } VecFunction createNext() { return new InfiniteChargedPlane(); } }; class InfiniteChargedPlane extends VecFunction { String getName() { return "infinite plane"; } void getField(double result[], double y[]) { double alpha = .0004; if (y[1] > -.01 && y[1] < .01) boundCheck = true; result[0] = 0; result[1] = (y[1] <= 0) ? alpha : -alpha; result[2] = java.lang.Math.abs(y[1])-1; } boolean checkBoundsWithForce() { return false; } VecFunction createNext() { return new Cylinder(); } }; class Cylinder extends VecFunction { String getName() { return "conducting cylinder"; } double getCylRadius() { return (aux1Bar.getValue()+1)/110.; } void setup() { setupBar(0, "Cylinder Size", 30); setupBar(1, "Cylinder Potential", 1); } void getField(double result[], double y[]) { // Smythe p70 double a = getCylRadius(); double farpt = 4; double pot = 2*(aux2Bar.getValue()/50.-1); final int mult = 4000; // add in charge that will bring cylinder to specified potential double cq = -pot/ (mult*(java.lang.Math.log(a)-java.lang.Math.log(farpt))); double pot0 = mult*cq*java.lang.Math.log(farpt); double y0 = y[0]; double y1 = y[1]; double r1 = distance(y0, y1); if (r1 < a) { result[0] = result[1] = 0; result[2] = pot; boundCheck = true; return; } double a1 = 5*cq/(r1*r1); result[0] = y0*a1; result[1] = y1*a1; result[2] = pot0-cq*mult*java.lang.Math.log(r1); } void calcDivergence() { double a = getCylRadius(); int i; for (i = 0; i != 100; i++) { double th = 2*pi*i/100.; double xx = java.lang.Math.cos(th)*a; double yy = java.lang.Math.sin(th)*a; grid[doubleToGrid(xx)][doubleToGrid(yy)].div -= reverse/20.; } } VecFunction createNext() { return new CylinderAndLineCharge(); } }; class CylinderAndLineCharge extends VecFunction { String getName() { return "cyl + line charge"; } double getCylRadius() { return (aux1Bar.getValue()+1)/110.; } double getSeparation() { return aux2Bar.getValue()/100.; } double getCylPos() { return getSeparation()/2; } double getPointPos() { return -getSeparation()/2-getCylRadius(); } void setup() { setupBar(0, "Cylinder Size", 30); setupBar(1, "Separation", 30); setupBar(2, "Cylinder Potential", 50); } void setupFrame() { // Smythe p70 q = -.0003; a = getCylRadius(); b = getSeparation() + a; spos = getCylPos(); imagePos = spos - a*a/b; double r2_0 = spos+a-imagePos; double r3_0 = spos+a-getPointPos(); // Find charge that will eliminate field at far end cylinder. // That's our definition of "ground". (We can't use the // potential at infinity as ground, since that diverges.) cq = a*a*(q/(r2_0*r2_0) - q/(r3_0*r3_0)); // calculate potential on cylinder with that charge. we will // subtract that value out of the potential everywhere. pot0 = -cq*java.lang.Math.log(a) +q*java.lang.Math.log(r2_0) -q*java.lang.Math.log(r3_0); // add in charge that will bring cylinder to specified potential cq -= (aux3Bar.getValue()/50.-1)*.0006/java.lang.Math.log(a); } double q, a, b, spos, imagePos, cq, pot0; void calcDivergence() { int i; double pos[] = rk_k1; double res[] = rk_k2; double a1 = getCylRadius()+.001; int x, y; for (x = 0; x != gridsize; x++) for (y = 0; y != gridsize; y++) grid[x][y].div = 0; grid[doubleToGrid(getPointPos())][doubleToGrid(0)].div = -reverse; for (i = 0; i != 200; i++) { double th = 2*pi*i/200.; double costh = java.lang.Math.cos(th); double sinth = java.lang.Math.sin(th); pos[0] = costh*a1+getCylPos(); pos[1] = sinth*a1; pos[2] = 0; curfunc.getField(res, pos); grid[doubleToGrid(costh*a+getCylPos())] [doubleToGrid(sinth*a)].div += (costh*res[0] + sinth*res[1])*60*reverse; } } void getField(double result[], double y[]) { // Smythe p70 double x1 = y[0]-spos; double r1 = distance(x1, y[1]); final int mult = 4000; double y0 = y[0]; double y1 = y[1]; if (r1 < a) { y0 = spos+a; y1 = 0; x1 = r1 = a; boundCheck = true; } double x2 = y0-imagePos; double r2 = distance(x2, y1); double x3 = y0-getPointPos(); double r3 = distance(x3, y1); double chargeSize = .001; if (r3 < chargeSize) boundCheck = true; double a1 = cq/(r1*r1); double a2 = -q/(r2*r2); double a3 = q/(r3*r3); result[0] = x1*a1 + x2*a2 + x3*a3; result[1] = y1*(a1+a2+a3); result[2] = mult*(-pot0 -cq*java.lang.Math.log(r1) +q*java.lang.Math.log(r2+1e-20) -q*java.lang.Math.log(r3+1e-20)); if (r1 == a) result[0] = result[1] = 0; } VecFunction createNext() { return new CylinderInField(); } }; class CylinderInField extends VecFunction { CylinderInField() { conducting = true; showD = false; } String getName() { return "cylinder in field"; } boolean conducting, showD; double a; void setupFrame() { a = aux1Bar.getValue()/100.; } void getField(double result[], double y[]) { // smythe p67 double a2 = a*a; double r = distance(y[0], y[1]); double e1 = aux2Bar.getValue()/10. + 1; double dimult = (conducting) ? 1 : (e1-1)/(e1+1); double fmult = .0006; final double potmult = 3; if (r < a) { result[0] = result[1] = result[2] = 0; if (conducting) boundCheck = true; else result[0] = (showD) ? e1*fmult*(1-dimult) : fmult*(1-dimult); result[2] = -potmult*(1-dimult)*y[0]; return; } double costh = y[0]/r; double sinth = y[1]/r; double r_2 = 1/(r*r); double er = (1+dimult*a2*r_2)*costh*fmult; double eth = -(1-dimult*a2*r_2)*sinth*fmult; er /= r; result[0] = y[0]*er - eth*sinth; result[1] = y[1]*er + eth*costh; result[2] = -potmult*(1-dimult*a2*r_2)*y[0]; } void setup() { setupBar(0, "Cylinder Size", 40); } void calcDivergence() { int i; double pos[] = rk_k1; double res[] = rk_k2; double a1 = a+.001; int x, y; for (x = 0; x != gridsize; x++) for (y = 0; y != gridsize; y++) grid[x][y].div = 0; for (i = 0; i != 200; i++) { double th = 2*pi*i/200.; double costh = java.lang.Math.cos(th); double sinth = java.lang.Math.sin(th); pos[0] = costh*a1; pos[1] = sinth*a1; pos[2] = 0; curfunc.getField(res, pos); double rx = res[0]; double ry = res[1]; double a2 = a-.001; pos[0] = costh*a2; pos[1] = sinth*a2; pos[2] = 0; curfunc.getField(res, pos); grid[doubleToGrid(costh*a)][doubleToGrid(sinth*a)].div += (costh*(rx-res[0]) + sinth*(ry-res[1]))*4e2*reverse; } } VecFunction createNext() { return new DielectricCylinderInFieldE(); } }; class DielectricCylinderInFieldE extends CylinderInField { DielectricCylinderInFieldE() { conducting = false; showD = false; } String getName() { return "dielec cyl in field"; } void setup() { setupBar(0, "Cylinder Size", 40); setupBar(1, "Dielectric Strength", 60); } VecFunction createNext() { return new SlottedPlane(); } }; /*class DielectricBoundaryE extends InverseSquaredRadial { boolean showD, conducting; double planeY; DielectricBoundaryE() { conducting = false; showD = false; } String getName() { return "dielec boundary"; } void setup() { setupBar(0, "Charge Location", 60); if (!conducting) setupBar(1, "Dielectric Strength", 60); } void getField(double result[], double y[]) { double cx = aux1Bar.getValue()/50.-1.001; double r1 = distance(y[0], y[1]-cx); if (r1 < chargeSize) boundCheck = true; double eps1 = 1; double eps2 = aux2Bar.getValue()/10. + 1; if (conducting) eps2 = 1e8; if (cx < planeY) { eps1 = eps2; eps2 = 1; } double q1 = .0003; double q2 = -(eps2-eps1)/(eps2+eps1)*q1; double ep = eps1; if (cx > planeY && y[1] < planeY || cx < planeY && y[1] > planeY) { q1 = 2*eps2*q1/(eps2+eps1); q2 = 0; ep = eps2; } double r2 = distance(y[0], y[1]-planeY*2+cx); result[2] = -2000*(q1/(r1*ep) + q2/(r2*ep)); if (!showD) { q1 /= ep; q2 /= ep; } double rq1 = q1/(r1*r1*r1); double rq2 = q2/(r2*r2*r2); result[0] = -y[0]*(rq1+rq2); result[1] = -(y[1]-cx)*rq1-(y[1]-planeY*2+cx)*rq2; } VecFunction createNext() { return new ConductingPlane(); } }; class ConductingPlane extends DielectricBoundaryE { ConductingPlane() { conducting = true; showD = false; planeY = -1; } String getName() { return "conducting plane + pt"; } VecFunction createNext() { return new SlottedPlane(); } };*/ class SlottedPlane extends VecFunction { String getName() { return "slotted conducting plane"; } Complex z, z2, z3; double getDivOffset() { return -17.3; } double getDivRange() { return 2.5; } SlottedPlane() { z = new Complex(); z2 = new Complex(); z3 = new Complex(); } void getField(double result[], double y[]) { // W = -.5E (z +- (z^2-a^2)^1/2) // dw/dz = -.5E (1 +- z/(z^2-a^2)^1/2) // Smythe p92 double a = (aux1Bar.getValue()+1)/101.; z.set(y[0], y[1]); if (y[1] >= -.01 && y[1] <= .01 && (y[0] < -a || y[0] > a)) { boundCheck = true; if (z.b == 0) z.b = -1e-8; } z2.set(z); z2.square(); z2.add(-a*a); z3.set(z2); z3.pow(.5); if (z3.b < 0) z3.mult(-1); z3.add(z.a, z.b); result[2] = z3.b*2; z2.pow(-.5); // I can barely understand what Smythe is talking about but // I think he wants the square root to always have a positive // imaginary part. Here we already took the reciprocal // (by calling pow(-.5) instead of pow(.5)) so we want the // root to have a negative imaginary part instead. if (z2.b > 0) z2.mult(-1); z2.mult(z); // field = (Im dw/dz, Re dw/dz) result[1] = -(1+z2.a) * .003; result[0] = -(z2.b) * .003; } void setup() { setupBar(0, "Slot Size", 30); } VecFunction createNext() { return new PlanePair(); } }; class PlanePair extends ConductingPlate { String getName() { return "conducting planes w/ gap"; } PlanePair() { super(); plate = false; } void setup() { setupBar(0, "Gap Size", 20); } double getDivOffset() { return -17; } VecFunction createNext() { return null; } }; #endif #ifndef BUILD_E class InverseRotational extends InverseRadial { String getName() { return BUILD_CASE_EMV(null, "current line", "1/r rotational"); } void getField(double result[], double y[]) { double r = distance(y[0], y[1]); if (showA) { result[0] = result[1] = 0; result[2] = -.001*(java.lang.Math.log(r)-.5); } else { if (r < lineWidth*2) boundCheck = true; rotateParticle(result, y, .0001/(r*r)); } } boolean nonGradient() { return true; } VecFunction createNext() { return new InverseRotationalPotential(); } }; class InverseRotationalPotential extends VecFunction { String getName() { return "1/r rotational potential"; } void getField(double result[], double y[]) { double r = distance(y[0], y[1]); rotateParticle(result, y, .0001/(r*r)); if (r < lineWidth*2) boundCheck = true; else if (y[0] >= 0 && y[1] < .001 && y[1] > -.025) { boundCheck = true; if (y[1] == 0) result[1] = 1e8; } double ang = java.lang.Math.atan2(y[1], y[0]); if (ang < 0) ang += 2*pi; result[2] = (pi-ang)*.3; } VecFunction createNext() { return new InverseRotationalDouble(); } }; class InverseRotationalDouble extends InverseRadialDouble { InverseRotationalDouble() { dir2 = 1; ext = false; } int dir2; boolean ext; String getName() { return BUILD_CASE_EMV(null, "current line double", "1/r rotational double"); } void getField(double result[], double y[]) { // we want to make sure sep is a multiple of the grid spacing // so that the divergence gets calculated correctly. double sep = gridToDouble(gridsize/2+ aux1Bar.getValue()*gridsize/200); double r = distance(y[0] - sep, y[1]); double r2 = distance(y[0] + sep, y[1]); if (ext) { double p = aux3Bar.getValue()*pi/50.; double s = aux2Bar.getValue()/6.; getDirectionField(result, y, p); result[0] *= s; result[1] *= s; result[2] *= s; } else result[0] = result[1] = result[2] = 0; if (showA) { if (dir2 == 1) result[2] += -.001*(java.lang.Math.log(r)+ java.lang.Math.log(r2)-1); else result[2] += .001*(java.lang.Math.log(r)- java.lang.Math.log(r2)); } else { if (r < lineWidth*2) boundCheck = true; rotateParticleAdd(result, y, .0001/(r*r), sep, 0); if (r2 < lineWidth*2) boundCheck = true; rotateParticleAdd(result, y, dir2*.0001/(r2*r2), -sep, 0); } } void setup() { setupBar(0, "Line Separation", 30); if (ext) { setupBar(1, "Ext. Strength", 7); setupBar(2, "Ext. Direction", 0); } } boolean nonGradient() { return true; } VecFunction createNext() { return new InverseRotationalDoubleExt(); } }; class InverseRotationalDoubleExt extends InverseRotationalDouble { InverseRotationalDoubleExt() { ext = true; } String getName() { return BUILD_CASE_EMV(null, "cur line double + ext", "1/r rot double + ext"); } VecFunction createNext() { return new InverseRotationalDipole(); } }; class InverseRotationalDipole extends InverseRotationalDouble { InverseRotationalDipole() { dir2 = -1; } String getName() { return BUILD_CASE_EMV(null, "current line dipole", "1/r rotational dipole"); } VecFunction createNext() { return new InverseRotationalDipoleExt(); } }; class InverseRotationalDipoleExt extends InverseRotationalDouble { InverseRotationalDipoleExt() { dir2 = -1; ext = true; } void setup() { super.setup(); aux2Bar.setValue(3); aux3Bar.setValue(25); } String getName() { return BUILD_CASE_EMV(null, "cur line dipole + ext", "1/r rot dipole + ext"); } VecFunction createNext() { return new OneDirectionFunction(); } }; class OneDirectionFunction extends VecFunction { String getName() { return BUILD_CASE_EMV(null, "uniform field", "one direction"); } void getField(double result[], double y[]) { double th = aux1Bar.getValue() * pi /50.; getDirectionField(result, y, th); } void setup() { setupBar(0, "Theta", 0); } VecFunction createNext() { return BUILD_CASE_EMV(null, new MovingChargeField(), new InverseSquaredRadialSphere()); } }; void getDirectionField(double result[], double y[], double th) { double sinth = java.lang.Math.sin(th); double costh = java.lang.Math.cos(th); if (!showA) { result[0] = .0003*costh; result[1] = .0003*sinth; result[2] = -.4*(y[0]*costh + y[1]*sinth); } else { // The A for this field is a linear rotation around // an axis pointing in the direction of the field. // First calculate the axis. double axis[] = new double[3]; axis[0] = costh; axis[1] = sinth; axis[2] = 0; // now project the point we are looking at onto the axis double d = dot(axis, y); // subtract out this projection so we can get the vector // from the point to the axis double r[] = new double[3]; int i; for (i = 0; i != 2; i++) r[i] = .0006*(y[i] - axis[i]*d); // cross this vector with the axis vector to get the result cross(result, axis, r); } } #endif #ifndef BUILD_V class MovingChargeField extends InverseSquaredRadial { String getName() { return "moving charge"; } void getField(double result[], double y[]) { double rz = distance(y); if (showA) { result[0] = result[1] = 0; result[2] = .0003/rz; } else { double r = distance(y[0], y[1]); if (rz < chargeSize) boundCheck = true; rotateParticle(result, y, .0001/(rz*rz*rz)); } } boolean nonGradient() { return true; } VecFunction createNext() { return BUILD_CASE_EMV(null, new FastChargeField(), null); } }; #endif #ifdef BUILD_V class InverseSquaredRadialSphere extends VecFunction { String getName() { return "1/r^2 sphere"; } double getSize() { return (aux1Bar.getValue()+1)/110.; } void getField(double result[], double y[]) { double r = distance(y); if (r < .01) boundCheck = true; double a = getSize(); result[2] = .2*((r > a) ? -1/r : -3/(2*a)+r*r/(2*a*a*a)); if (r < a) r = a; double alpha = .0003/(r*r*r); result[0] = -y[0]*alpha; result[1] = -y[1]*alpha; } void setup() { setupBar(0, "Sphere Size", 30); } boolean checkBoundsWithForce() { return false; } VecFunction createNext() { return new ConstRadial(); } }; class ConstRadial extends InverseSquaredRadial { String getName() { return "const radial"; } void getField(double result[], double y[]) { double r = distance(y); if (r < chargeSize) boundCheck = true; double q = .0006/r; result[0] = -q*y[0]; result[1] = -q*y[1]; result[2] = r; } boolean checkBoundsWithForce() { return false; } VecFunction createNext() { return new LinearRadial(); } }; class LinearRadial extends InverseSquaredRadial { String getName() { return "linear radial"; } void getField(double result[], double y[]) { double r = distance(y); if (r < chargeSize) boundCheck = true; double k = .0009; result[0] = -y[0]*k; result[1] = -y[1]*k; result[2] = r*r-1; } boolean checkBoundsWithForce() { return false; } VecFunction createNext() { return new ConstantToYAxis(); } }; class ConstantToYAxis extends InverseRadial { String getName() { return "constant to y axis"; } void getField(double result[], double y[]) { double alpha = .0004; if (y[0] > -.01 && y[0] < .01) boundCheck = true; result[0] = (y[0] <= 0) ? alpha : -alpha; result[1] = 0; result[2] = java.lang.Math.abs(y[0])-1; } boolean checkBoundsWithForce() { return false; } VecFunction createNext() { return new LinearToYAxis(); } }; class LinearToYAxis extends InverseRadial { String getName() { return "linear to y axis"; } void getField(double result[], double y[]) { double r = java.lang.Math.abs(y[0]); if (r < lineWidth) boundCheck = true; double q = .0009; result[0] = -y[0]*q; result[1] = 0; result[2] = r*r-1; } boolean checkBoundsWithForce() { return false; } VecFunction createNext() { return new LinearToXYAxes(); } }; class LinearToXYAxes extends VecFunction { String getName() { return "2-dimensional oscillator"; } void getField(double result[], double y[]) { double alpha = .0006; double r = java.lang.Math.sqrt((aux1Bar.getValue()+1)/51.); result[0] = -alpha*r*y[0]; result[1] = -alpha/r*y[1]; result[2] = (y[0]*y[0]*r+y[1]*y[1]/r)-1; } void setup() { setupBar(0, "X/Y Ratio", 50); } boolean checkBoundsWithForce() { return false; } VecFunction createNext() { return new InverseToYAxis(); } }; class InverseToYAxis extends VecFunction { String getName() { return "inverse to y axis"; } void getField(double result[], double y[]) { if (y[0] > -.01 && y[0] < .01) boundCheck = true; double alpha = .0003; double zz = y[0]; if (zz == 0) zz = .00001; result[0] = -alpha/zz; result[1] = 0; result[2] = -.01/(zz*zz); } VecFunction createNext() { return new InverseSquareRotational(); } }; class InverseSquareRotational extends InverseRadial { String getName() { return "1/r^2 rotational"; } void getField(double result[], double y[]) { double r = distance(y[0], y[1]); if (r < lineWidth*2) boundCheck = true; rotateParticle(result, y, .0001/(r*r*r)); } boolean nonGradient() { return true; } VecFunction createNext() { return new LinearRotational(); } }; class LinearRotational extends InverseRadial { String getName() { return "linear rotational"; } void getField(double result[], double y[]) { double q = .0006; result[0] = -q*y[1]; result[1] = q*y[0]; result[2] = 0; } boolean nonGradient() { return true; } VecFunction createNext() { return new ConstantRotational(); } }; class ConstantRotational extends InverseRadial { String getName() { return "constant rotational"; } void getField(double result[], double y[]) { double r = distance(y[0], y[1]); rotateParticle(result, y, .0006/r); } boolean nonGradient() { return true; } VecFunction createNext() { return new FxEqualsYField(); } }; class FxEqualsYField extends VecFunction { String getName() { return "(y,0)"; } void getField(double result[], double y[]) { result[2] = result[1] = 0; result[0] = y[1] * .0009; } boolean nonGradient() { return true; } VecFunction createNext() { return new FxEqualsY2(); } }; class FxEqualsY2 extends VecFunction { String getName() { return "(y^2,0)"; } void getField(double result[], double y[]) { result[2] = result[1] = 0; result[0] = y[1]*y[1] * .002; } boolean nonGradient() { return true; } VecFunction createNext() { return new Saddle(); } }; class Saddle extends VecFunction { String getName() { return "saddle"; } void getField(double result[], double y[]) { double q = .001; result[0] = -q*y[0]; result[1] = q*y[1]*.5; result[2] = 1*(2*y[0]*y[0] - y[1]*y[1]); } VecFunction createNext() { return new RotationalExpansion(); } }; class RotationalExpansion extends VecFunction { String getName() { return "rotation + expansion"; } void getField(double result[], double y[]) { double q = .0006; result[0] = q*(y[0]-y[1]); result[1] = q*(y[0]+y[1]); result[2] = 0; } boolean nonGradient() { return true; } VecFunction createNext() { return new Function4Field(); } }; class Function4Field extends VecFunction { String getName() { return "(x^2-y,x+y^2)"; } void getField(double result[], double y[]) { double q = .0003; result[0] = q*(y[0]*y[0]-y[1]); result[1] = q*(y[0]+y[1]*y[1]); result[2] = 0; } boolean nonGradient() { return true; } VecFunction createNext() { return new Function5Field(); } }; class Function5Field extends VecFunction { String getName() { return "(x+y^2,x^2-y)"; } void getField(double result[], double y[]) { double q = .0003; result[0] = q*(y[0]+y[1]*y[1]); result[1] = q*(y[0]*y[0]-y[1]); result[2] = 0; } boolean nonGradient() { return true; } VecFunction createNext() { return new Function6Field(); } }; class Function6Field extends VecFunction { String getName() { return "(x,x^2)"; } void getField(double result[], double y[]) { double q = .0006; result[0] = q*y[0]; result[1] = q*y[0]*y[0]; result[2] = 0; } boolean nonGradient() { return true; } VecFunction createNext() { return new Function7Field(); } }; class Function7Field extends VecFunction { String getName() { return "u=x^2+y"; } void getField(double result[], double y[]) { double q = .0003; result[0] = -2*y[0]*q; result[1] = -q; result[2] = (y[0]*y[0]+y[1])*1; } VecFunction createNext() { return new PendulumPotential(); } }; class PendulumPotential extends VecFunction { String getName() { return "pendulum potential"; } void getField(double result[], double y[]) { double q = .0006; double xx = y[0]*3.1; double yy = y[1]*3.1; double cosx = java.lang.Math.cos(xx); double cosy = java.lang.Math.cos(yy); double sinx = java.lang.Math.sin(xx); double siny = java.lang.Math.sin(yy); result[0] = -q*sinx*cosy; result[1] = -q*cosx*siny; result[2] = -cosx*cosy*.5; } VecFunction createNext() { return new Function8Field(); } }; class Function8Field extends VecFunction { String getName() { return "sin(r2)/r2"; } void getField(double result[], double y[]) { double r2 = 23*(y[0]*y[0]+y[1]*y[1])+.00000001; result[2] = java.lang.Math.sin(r2)/r2; double r= java.lang.Math.sqrt(r2); double sinr2 = java.lang.Math.sin(r2); double cosr2 = java.lang.Math.cos(r2); double r4 = r2*r2; double vecR = (sinr2/r4-cosr2/r2)*.01; result[0] = y[0]*vecR; result[1] = y[1]*vecR; } VecFunction createNext() { return new UserDefinedPotential(); } }; class UserDefinedPotential extends VecFunction { Expr expr; double y0[]; String getName() { return "user-defined potential"; } void setup() { textFields[0].setText("x*x"); textFields[0].show(); textFieldLabel.setText("Potential Function"); textFieldLabel.show(); actionPerformed(); y0 = new double[3]; } void actionPerformed() { parseError = false; ExprParser ep = new ExprParser(textFields[0].getText()); expr = ep.parseExpression(); if (ep.gotError()) parseError = true; functionChanged = true; } void getField(double result[], double y[]) { double k = .00001; int i; for (i = 0; i != 3; i++) y0[i] = y[i]; double pot0 = expr.eval(y0); y0[0] += k; result[0] = pot0-expr.eval(y0); y0[0] = y[0]; y0[1] += k; result[1] = pot0-expr.eval(y0); y0[1] = y[1]; result[2] = pot0*.01; for (i = 0; i != 2; i++) if (!(result[i] > -10 && result[i] < 10)) boundCheck = true; } VecFunction createNext() { return new UserDefinedFunction(); } }; class UserDefinedFunction extends VecFunction { Expr exprs[]; String getName() { return "user-defined field"; } void setup() { exprs = new Expr[3]; textFields[0].setText("x"); textFields[1].setText("y"); textFieldLabel.setText("Field Functions"); textFieldLabel.show(); int i; for (i = 0; i != 2; i++) textFields[i].show(); actionPerformed(); } void actionPerformed() { int i; parseError = false; for (i = 0; i != 2; i++) { ExprParser ep = new ExprParser(textFields[i].getText()); exprs[i] = ep.parseExpression(); if (ep.gotError()) parseError = true; } functionChanged = true; } void getField(double result[], double y[]) { double k = .0002; int i; for (i = 0; i != 2; i++) { result[i] = k*exprs[i].eval(y); if (!(result[i] > -10 && result[i] < 10)) boundCheck = true; } result[2] = 0; } boolean nonGradient() { return true; } VecFunction createNext() { return null; } }; #endif class DrawData { public Graphics g; public double mult; public double field[], vv[]; }; class Particle { public double pos[]; public double vel[]; public double lifetime; public double phi, theta, phiv, thetav; public double stepsize; public Color color; Particle() { pos = new double[3]; vel = new double[3]; stepsize = 1; } }; class FieldVector { public int sx1, sy1, sx2, sy2; public double p1[], p2[]; public int col; public int viewPri; }; class GridElement { public double height, div, curl, normdot, vecX, vecY; public boolean visible, valid; }; #ifdef BUILD_E // this is a complex number class that avoids creating new objects // whenever possible. Takes a little more effort to use, but it is // much faster because we don't have to do constant garbage // collection. class Complex { public double a, b; Complex() { a = b = 0; } void set(double aa, double bb) { a = aa; b = bb; } void set(Complex c) { set(c.a, c.b); } void add(double r) { a += r; } void add(double r, double i) { a += r; b += i; } void square() { set(a*a-b*b, 2*a*b); } void mult(double c, double d) { set(a*c-b*d, a*d+b*c); } void mult(double c) { a *= c; b *= c; } void mult(Complex c) { mult(c.a, c.b); } void recip() { double n = a*a+b*b; set(a/n, -b/n); } void pow(double p) { double arg = java.lang.Math.atan2(b, a); arg *= p; double abs = java.lang.Math.pow(a*a+b*b, p*.5); set(abs*java.lang.Math.cos(arg), abs*java.lang.Math.sin(arg)); } void sin() { set(cosh(b)*java.lang.Math.sin(a), java.lang.Math.cos(a)*sinh(b)); } void cos() { set(cosh(b)*java.lang.Math.cos(a), java.lang.Math.sin(a)*sinh(b)); } void log() { set(java.lang.Math.log(a*a+b*b), java.lang.Math.atan2(b, a)); } void arcsin() { Complex z2 = new Complex(); z2.set(a, b); z2.square(); z2.mult(-1); z2.add(1); z2.pow(.5); mult(0, 1); add(z2.a, z2.b); log(); mult(0, -1); } }; double cosh(double a) { return .5*(java.lang.Math.exp(a)+java.lang.Math.exp(-a)); } double sinh(double a) { return .5*(java.lang.Math.exp(a)-java.lang.Math.exp(-a)); } #endif }; #ifdef BUILD_V class ExprState { public double x, y, z; } class Expr { Expr(Expr e1, Expr e2, int v) { left = e1; right = e2; type = v; } Expr(int v, double vv) { type = v; value = vv; } Expr(int v) { type = v; } double eval(double es[]) { switch (type) { case E_ADD: return left.eval(es)+right.eval(es); case E_SUB: return left.eval(es)-right.eval(es); case E_MUL: return left.eval(es)*right.eval(es); case E_DIV: return left.eval(es)/right.eval(es); case E_POW: return java.lang.Math.pow(left.eval(es), right.eval(es)); case E_UMINUS: return -left.eval(es); case E_VAL: return value; case E_X: return es[0]*10; case E_Y: return es[1]*10; case E_R: return java.lang.Math.sqrt( es[0]*es[0]+es[1]*es[1])*10; case E_SIN: return java.lang.Math.sin(left.eval(es)); case E_COS: return java.lang.Math.cos(left.eval(es)); case E_ABS: return java.lang.Math.abs(left.eval(es)); case E_EXP: return java.lang.Math.exp(left.eval(es)); case E_LOG: return java.lang.Math.log(left.eval(es)); case E_SQRT: return java.lang.Math.sqrt(left.eval(es)); case E_TAN: return java.lang.Math.tan(left.eval(es)); default: System.out.print("unknown\n"); } return 0; } Expr left, right; double value; int type; static final int E_ADD = 1; static final int E_SUB = 2; static final int E_X = 3; static final int E_Y = 4; static final int E_Z = 5; static final int E_VAL = 6; static final int E_MUL = 7; static final int E_DIV = 8; static final int E_POW = 9; static final int E_UMINUS = 10; static final int E_SIN = 11; static final int E_COS = 12; static final int E_ABS = 13; static final int E_EXP = 14; static final int E_LOG = 15; static final int E_SQRT = 16; static final int E_TAN = 17; static final int E_R = 18; }; class ExprParser { String text; String token; int pos; int tlen; boolean err; void getToken() { while (pos < tlen && text.charAt(pos) == ' ') pos++; if (pos == tlen) { token = ""; return; } int i = pos; int c = text.charAt(i); if ((c >= '0' && c <= '9') || c == '.') { for (i = pos; i != tlen; i++) { if (!((text.charAt(i) >= '0' && text.charAt(i) <= '9') || text.charAt(i) == '.')) break; } } else if (c >= 'a' && c <= 'z') { for (i = pos; i != tlen; i++) { if (!(text.charAt(i) >= 'a' && text.charAt(i) <= 'z')) break; } } else { i++; } token = text.substring(pos, i); pos = i; } boolean skip(String s) { if (token.compareTo(s) != 0) return false; getToken(); return true; } void skipOrError(String s) { if (!skip(s)) err = true; } Expr parseExpression() { if (token.length() == 0) return new Expr(Expr.E_VAL, 0.); Expr e = parse(); if (token.length() > 0) err = true; return e; } Expr parse() { Expr e = parseMult(); while (true) { if (skip("+")) e = new Expr(e, parseMult(), Expr.E_ADD); else if (skip("-")) e = new Expr(e, parseMult(), Expr.E_SUB); else break; } return e; } Expr parseMult() { Expr e = parseUminus(); while (true) { if (skip("*")) e = new Expr(e, parseUminus(), Expr.E_MUL); else if (skip("/")) e = new Expr(e, parseUminus(), Expr.E_DIV); else break; } return e; } Expr parseUminus() { skip("+"); if (skip("-")) return new Expr(parsePow(), null, Expr.E_UMINUS); return parsePow(); } Expr parsePow() { Expr e = parseTerm(); while (true) { if (skip("^")) e = new Expr(e, parseTerm(), Expr.E_POW); else break; } return e; } Expr parseFunc(int t) { skipOrError("("); Expr e = parse(); skipOrError(")"); return new Expr(e, null, t); } Expr parseTerm() { if (skip("(")) { Expr e = parse(); skipOrError(")"); return e; } if (skip("x")) return new Expr(Expr.E_X); if (skip("y")) return new Expr(Expr.E_Y); if (skip("r")) return new Expr(Expr.E_R); if (skip("pi")) return new Expr(Expr.E_VAL, 3.14159265358979323846); if (skip("e")) return new Expr(Expr.E_VAL, 2.7182818284590452354); if (skip("sin")) return parseFunc(Expr.E_SIN); if (skip("cos")) return parseFunc(Expr.E_COS); if (skip("abs")) return parseFunc(Expr.E_ABS); if (skip("exp")) return parseFunc(Expr.E_EXP); if (skip("log")) return parseFunc(Expr.E_LOG); if (skip("sqrt")) return parseFunc(Expr.E_SQRT); if (skip("tan")) return parseFunc(Expr.E_TAN); try { Expr e = new Expr(Expr.E_VAL, Double.valueOf(token).doubleValue()); getToken(); return e; } catch (Exception e) { err = true; System.out.print("unrecognized token: " + token + "\n"); return new Expr(Expr.E_VAL, 0); } } ExprParser(String s) { text = s; tlen = text.length(); pos = 0; err = false; getToken(); } boolean gotError() { return err; } }; #endif interface DecentScrollbarListener { abstract void scrollbarValueChanged(DecentScrollbar ds); abstract void scrollbarFinished(DecentScrollbar dc); } // this is a scrollbar that notifies us when the user is _done_ fiddling // with the value. class DecentScrollbar extends Canvas implements MouseListener, MouseMotionListener { int value, lo, hi; DecentScrollbarListener listener; DecentScrollbar(DecentScrollbarListener parent, int start, int lo_, int hi_) { value = start; lo = lo_; hi = hi_; listener = parent; addMouseListener(this); addMouseMotionListener(this); gray1 = new Color(104, 104, 104); gray2 = new Color(168, 168, 168); gray3 = new Color(192, 192, 192); gray4 = new Color(224, 224, 224); } Color gray1, gray2, gray3, gray4; public Dimension getPreferredSize() { return new Dimension(20,20); } static final int tw = 8; boolean dragging; int thumbpos, dragoffset; public void paint(Graphics g) { Dimension size = getSize(); int w = size.width; int h = size.height; int x = thumbpos = (value-lo)*(w-2-tw)/(hi-lo)+1; g.setColor(gray2); g.fillRect(0, 0, w, h); g.setColor(gray3); g.fillRect(x, 2, tw, h-4); g.setColor(gray4); g.drawLine(0, h-1, w, h-1); g.drawLine(w-1, 1, w-1, h-1); g.drawLine(x, 1, x+tw-1, 1); g.drawLine(x, 1, x, h-2); g.setColor(gray1); g.drawLine(0, 0, w-1, 0); g.drawLine(0, 0, 0, h-1); g.drawLine(x+tw-1, 2, x+tw-1, h-2); g.drawLine(x+1, h-2, x+tw-1, h-2); } int getValue() { return value; } boolean setValue(int v) { if (v < lo) v = lo; if (v > hi) v = hi; if (value == v) return false; value = v; repaint(); return true; } public void mousePressed(MouseEvent e) { if (thumbpos <= e.getX() && thumbpos+tw >= e.getX()) { dragging = true; dragoffset = e.getX()-thumbpos; } } public void mouseReleased(MouseEvent e) { if (dragging) listener.scrollbarFinished(this); dragging = false; } public void mouseClicked(MouseEvent e) { } public void mouseEntered(MouseEvent e) { } public void mouseExited(MouseEvent e) { } public void mouseDragged(MouseEvent e) { if (!dragging) return; int x = e.getX()-dragoffset; Dimension size = getSize(); int v = (x-1)*(hi-lo)/(size.width-2-tw)+lo; if (setValue(v)) listener.scrollbarValueChanged(this); } public void mouseMoved(MouseEvent e) { } }; A far higher place must be assigned to Judaism among the competitors for the allegiance of Europe. The cosmopolitan importance at one time assumed by this religion has been considerably obscured, owing to the subsequent devolution of its part to Christianity. It is, however, by no means impossible that, but for the diversion created by the Gospel, and the disastrous consequences of their revolt against Rome, the Jews might have won the world to a purified form of their own monotheism. A few significant circumstances are recorded showing how much influence they had acquired, even in Rome, before the first preaching of Christianity. The first of these is to be found in Cicero¡¯s defence of Flaccus. The latter was accused of appropriating part of the annual contributions sent to the temple at Jerusalem; and, in dealing with this charge, Cicero speaks of the Jews, who were naturally prejudiced against his client, as a powerful faction the hostility of which he is anxious not to provoke.330 Some twenty years later, a great advance has been made. Not only must the material interests of the Jews be respected, but a certain conformity to their religious prescriptions is considered a mark of good breeding, In one of his most amusing satires, Horace tells us how, being anxious to shake off a bore, he appeals for help to his friend Aristius Fuscus, and reminds him of217 some private business which they had to discuss together. Fuscus sees his object, and being mischievously determined to defeat it, answers: ¡®Yes, I remember perfectly, but we must wait for some better opportunity; this is the thirtieth Sabbath, do you wish to insult the circumcised Jews?¡¯ ¡®I have no scruples on that point,¡® replies the impatient poet. ¡®But I have,¡¯ rejoins Fuscus,¡ª¡®a little weak-minded, one of the many, you know¡ªexcuse me, another time.¡®331 Nor were the Jews content with the countenance thus freely accorded them. The same poet elsewhere intimates that whenever they found themselves in a majority, they took advantage of their superior strength to make proselytes by force.¡¯332 And they pursued the good work to such purpose that a couple of generations later we find Seneca bitterly complaining that the vanquished had given laws to the victors, and that the customs of this abominable race were established over the whole earth.333 Evidence to the same effect is given by Philo Judaeus and Josephus, who inform us that the Jewish laws and customs were admired, imitated, and obeyed over the whole earth.334 Such assertions might be suspected of exaggeration, were they not, to a certain extent, confirmed by the references already quoted, to which others of the same kind may be added from later writers showing that it was a common practice among the Romans to abstain from work on the Sabbath, and even to celebrate it by praying, fasting, and lighting lamps, to visit the synagogues, to study the law of Moses, and to pay the yearly contribution of two drachmas to the temple at Jerusalem.335 Jeff¡¯s hand was quietly coming down. ¡°What happened to you?¡± begged Sandy. ¡°Something new has come up, sir. I was waiting there by my ship a good while back, and I heard another one cruising and spiraling, shooting the field, I guess, because he came in and set down. My crate, just the way you ordered, was down by the grove, not in plain sight in the middle of the course. But Jeff set his ship down, left the engine running, and went off. I stayed hid to see what would happen, but when he didn¡¯t come back, I thought I¡¯d better go and find you¡ªand see if it meant anything to you.¡± She stood alone, with the sticky, wet knife in her hand, catching her breath, coming out of the madness. Then she stooped, and pushing the branches aside felt about for her pistol. It lay at the root of a tree, and[Pg 80] when she had picked it up and put it back in the holster, there occurred to her for the first time the thought that the shot in the dead stillness must have roused the camp. And now she was sincerely frightened. If she were found here, it would be more than disagreeable for Landor. They must not find her. She started at a swift, long-limbed run, making a wide detour, to avoid the sentries, bending low, and flying silently among the bushes and across the shadowy sands. The year 1756 opened with menaces to England of the most serious nature. The imbecility of the Ministry was beginning to tell in the neglect of its colonies and its defences. France threatened to invade us, and a navy of fifty thousand men was suddenly voted, and an army of thirty-four thousand two hundred and sixty-three of native troops; but as these were not ready, it was agreed to bring over eight thousand Hessians and Hanoverians. To pay for all this it was necessary to grant excessive supplies, and lay on new duties and taxes. In presenting the money bills in the month of May, Speaker Onslow could not avoid remarking that there were two circumstances which tended to create alarm¡ªforeign subsidies and foreign troops introduced, and nothing but their confidence in his Majesty could allay their fears, or give them confidence that their burdens would be soon reduced. There was, in fact, no chance for any such reduction, for wars, troubles, and disgraces were gathering around from various quarters. The first reverse came from the Mediterranean. MUCH to their amazement, the boys waked up the next morning in Nashville, and found that they had passed through the "dark and bloody ground" of Kentucky absolutely without adventure. After drawing and dividing the rations and cartridges. Si gave the boys the necessary instruction about having their things ready so that they could get them in the dark the next morning, and ordered them to disregard the bonfires and mirth-making, and lie down to get all the sleep they could, in preparation for the hard work of the next day. Then, like the rest of the experienced men, who saw that the campaign was at length really on, and this would be the last opportunity for an indefinite while to write, he sat down to write short letters to his mother and to Annabel. "Bully for the Wild Wanderers of the Wabash," Shorty joined in. "They're the boss regiment in the army o' the Cumberland, and the Army o' the Cumberland's the boss army on earth. Hooray for US Co. Le's have a speech. Where's Monty Scruggs?" "Bring a light, do¡ªI can't abide this dark." Albert suddenly began to look uneasy. After all he was not really drunk, only a little fuddled. He walked straight, and his roll was natural to him, while though he was exceedingly cheerful, and often burst into song, his words were not jumbled, and he generally seemed to have a fair idea of what he was saying. "But I heard what the doctor said to you." "A purty accident¡ªwud them stacks no more dry than a ditch. 'Twas a clear case of 'bustion¡ªfireman said so to me; as wicked and tedious a bit o' wark as ever I met in my life." Calverley stept from the shadow of the cliff, and beheld a meteor in the sky, brightening and expanding, as the clouds opened, until it assumed the appearance of a brilliant star, of astonishing magnitude, encircled by dazzling rays, which, in a singular manner, were all inclined in one direction, and pointing to that part of the horizon where lay the rival of England¡ªFrance. The foreman's face assumed a deeper hue than usual: he looked fiercely at the galleyman, but there was a determination in the weather-beaten face that made him pause ere he spoke. "Galleyman," he at length said, "you knew the business before you came: if you be so fond of saving old witches' lives, why didn't you say so, that I might not now be in this dilemma?" "No, no, not the boy," replied Merritt, rather impatiently. HoMEÁíÀàС˵ߣɫ ENTER NUMBET 0017
pddchain.com.cn
tehe9.com.cn
www.sakka.com.cn
xl6666.com.cn
www.daren2.net.cn
jiuqu6.net.cn
lichu5.com.cn
75ze.com.cn
09907.com.cn
1genzhen.com.cn
美国色迷迷导航 和美女打泡经历 里诺明日香 zui淫当的女人 欧美露逼网 WWW_233JJ_COM 色色爽图 裸体婚沙照图片 WWW.787ZZ.COM WWW.LLDM.COM WWW.CCC403.COM WWW.T262.COM WWW.WDHGB.COM WWW.YY9PLAY.COM WWW.NNN89.COM WWW.RIJIALU6.COM WWW.99DAO.COM WWW.44NANA.COM WWW.XXXGAO.COM WWW.CCC981.COM WWW.JINGHAI.COM WWW.1024GC.COM WWW.999338.COM WWW.PV990.COM WWW.BANGGOU.COM WWW.ANQIMA.COM WWW.UZ92.COM LL.COOL WWW.JTBZJZ.NET WWW.SHJRCS.COM WWW.S5CM.COM WWW.WXKPSS.COM WWW.G5RZ.COM WWW.WH-YDNT.COM 出轨的乱伦小说 欧美性爱喷 lubi9999 先锋熟女人妻第一页 亚州依依影院 黄色慰安妇电影 校园激情卡通贴图 操性感人妻照片 90后色 免费视频自拍做爱 在线播放bbb39卡顿无需播放器 公马母毛驴交配 意淫强奸人妻女友 手机在线免费av梦工厂a片 巨乳肥臀15P 迷奸漫画 wwwtiantang 色图777 色网姐姐五月 av12-av在线日本视频在线观看 偷拍自拍做爱小说 看大片wwwbbb330com videosgratatistv欲女 巴西美女三级性爱 西瓜操 97超碰乱伦强奸资源视 国产自拍91gav 免费在线做爱网站ed2k 幼女啪啪啪影院 mmm9ckanzycom 强奸乱伦自慰丝袜亚洲色图欧美电影偷拍自拍在线电影好友自慰热の 狂威杀手 黄漫视频网址黄漫视频网址大全 啊啊啊骚穴好痒插 西西人艺体艺摄影 免费成人在线视频导航 短篇合集 chinesepornvideos--porn300 wwwya778com 有码伦理 8090狠狠射 香港三级片台三级片港台三级影院 性吧春天 中国情人黄黄色网站 wwwaaaa49n 三汲片电影视频播放 色姐妹与狗 影音先锋资源av成人片 aaa黄色日本网站播 亚洲色图欧美色图丰乳肥臀色图 77pian_com 英文字典 超碰在线看视频 嫩汁50p 高清av影院点击av 一淫妻 包射丝袜 偷拍自拍清纯唯美先锋 www成人综合网站com 清纯唯美88ttile88 妹妹给我下安眠药全节 91pro内部地址 淫姐姐影院乱伦小说 顶级裸体女模特 12jiusetengcom WWW756SqCOm JB123456 色色视频下载 作者不详bt工厂 日本乱伦AV下载下载 俄罗斯姑娘阴毛 岛国AV无码 色久妹在线视频 老和尚成人影院 啊嗯用力嗯好大动态图 幼女雄鲍 受法律堡海外mp4 幼h视频 艳母无需插件播放的网站 av作品丝袜恋足番号 家挺淫乱系列 姑姑的 哪里看av影片 呜呜呜呜色色色 操韩国空姐 猛男以超快速度不间断最后干的女友高潮长叫一声 做爱动态图片黄色片 田野麻衣 黄色片哥哥妹妹 玲花裸体照 97sec1 与英语老师做爱的小说 惠美里大战黑人 滛漫画 亚洲10p小穴 很很撸乱伦小说 黑人的大鸡巴大战波多野 李宗瑞mov275g 李宗瑞qq群论坛 taosenvcom 中年裸乳 欧美色图 迅雷下载 媚娘三洞被插 乳体兆片 父和女儿肏屄乱伦小说 欧美7岁幼女性交小电影 我爱x逼网 丝袜制服插插色色色 xxoo合成照 能播放的欧美群交视频 中文网成人电影 百度影音能打开的黄网 性感邻居 电影 qvod女无码 成人午夜快播百度影音 夫妻黄色三级片 我姐尻屁片 玩北京姐姐水穴 你们操逼喜欢摸奶吗操操 少妇露阴西西人体 爆操60岁老女人路脸视频 湖南省人民检察院 哪种网站你懂的 性吧欧美色图 动画版裸体美女图片 公公尻逼图片 强奸巨乳妇 筱声色色 yqazhouseqing 越南大胆室内人体艺术 黄图paoniu 亚洲日本色露逼 人体阴部裸图 你jj是否也硬了起来【24p】 国产自拍图片区国产区小说区 shaofunen射15p 妈妈被我一尺多长的大鸡巴干昏了 撸熟女骚 美女双飞国语 WWW187KKKCOM 亚寂寞美女用手指抠逼草莓 亚洲色图少妇做美图 黄色小说全集百度云 少妇流精视频 插屁眼av 领久美女热图 蜜穴偷情删除 WWWID8888COMYKT 美女如玉人体艺术 奥州性爱电影 哪里有幼女裸图网站 品se堂永久免费的论坛 干嫂嫂网站 苍井空艳照写真 富婆与动物乱伦 zaixianshipintaiwan www操逼聪明 做爱玩15p 少妇骚粉红色屄图 日本天天草莓影音 超碰白丝 肥婆肛门交 兽幼怀孕 日本儿子趁妈妈洗澡时强奸妈妈 草莓100加多撸 波多野结衣老师动图 脱衣女奴 印度毛片ed2k 中国大妈自拍 37丝袜资源论坛 色狼撸管网 久久精品视频在线看成人动漫 q2002影视院magnet 蕙和惠 qv0d快播电影妹妹嫩屄 国外幼幼资源站 日本少归内衣 看免黄色电影和黄色小说 美女大胆掰开阴道图片 小穴好嫩电影在线 wwwwhaorenshuocom avHS动漫网站 何洁色 9先锋资源com 小泽玛利亚高潮图片搜索 色狗成人激情网 逍遥坊成人电影 hentiahaven AV色图网站 美女bb插jj上流水 畅想听吧有声小说 红花樱井莉亚 小泽玛利亚天空 没封的h网 电脑加www打 开心妞妞五月天 怎么在手机上看黄片 如何看黄片 日韩黄色小说 色贝贝电影 小色孩导航 美女高清电影 日日啪 撸踏踏 撸飞哥 911色主站 voyeur浴场 fdd 1231av女主演 av页面导航 av精品 电车痴汉 新手村伦理影视 大哥色导航 - 百度 全裸激情插逼视屏 汤先生约炮童颜巨乳美女近景拍摄2 我的二婶韩国伦理电影 青青私密视频 人狗资源福利在线 日本校花黑丝视频 日本一本道高清无码a 日本男生插日本女生屁股视频 日韩口暴吞精在线视频 影音先锋观看仁科百华电影 淫穴搜 日本推油无码av magnet 日本性感诱惑少妇视频 日本最变态的绑法视频 自拍快播 妻子美脚连裤在线观看 青少年阴蒂毛茸茸视频 张柏芝32分视频神马影院 91视频播放列表 冲田杏梨wanz265在线看 selaosiji 台湾一级毛片下载 日韩伦理片制服丝袜 黄色视频美女丝袜 我一拍福利视频 欧美大屁股迅雷下载 ed2k 大香蕉一本道富二代 母亲和儿子做爱的伦理片 另类激情小说 av按摩 欧美足球av番号迅雷下载链接 www4445f 麻酥酥哟完整视频 magnet 国产草莓97在线 安卓手机免费在线黄片 爆乳i罩杯女教师高速夹射在线播放 撸管操逼影院 色爱亚洲五月天 调教青楼虐阴小说 AV天堂网 万能艺图网 精品av集锦 美腿紫竹领 半曰轻狂记手抄本 足控电影有哪些 扒衣云中文字幕 小米手工粉加盟jv168 (女装风衣) -(男休闲皮鞋) gvg一173在线观看 BB成人影院 瑞逼综合成人网 草民影院偷拍自拍无码 xxo影院 XXOO84短视频 torrenkity影音先锋 口工日本无翼岛视频全集 猎艳自拍 国产经典 日本操逼带声音 波儿影院 259LUXU-920 mp4 雾岛樱自慰 迅雷无码张柏芝种子 东方久久热 年轻的岳毌2 偷拍自拍爽图 茱莉亚·文斯H 91色老板福利电影 在线黄片主播。 一本道京东热手机视频在线 最新日韩性爱视频 www98bncom 天天干日夜谢 草莓论坛 se 吻戏软件影院 裸五月秋色 正在播放国产优衣库 4438xxx16 一本道av无码1v1爱泽花梨 李宗瑞在线偷拍视频 泷泽萝拉在线云盘 冲田杏梨爱剪辑在线无码观看 高手解玄区一好运一点通 操逼人兽com 周末同床 大尺度高清无码视频网站 1邪恶帝国漫画 日本女同 影音先锋男人资源站 日本偷拍厕所小便视频 迅雷下载 AV酒井法子黄片 亚美成人三级片 俺也去俺也干色啪啪 北川童在线 屌麻B视频 拍拍影院三级片 于子涵福利视频在线 小夫妻和单男玩剧情扮演 限制片zaixianbofang 小仙儿叫床声 性爱A V视频 国产三级浴室 桃色影院爱啪网 日本东凛视频在线 噜噜网噜噜巴噜噜色午夜影院 汤姆影院tam四虎影库 坐爱视一 李雅在线艺泽 淫乱丝袜 老师和学生腿上挂了安全套里番 上榴社区续缘舍 自拍偷拍日韩视频国产自拍 台湾妹午夜福利在线观看 汤姆猫成人影院的网址 谷露三级片 东京强制系列免费观看 操逼视频啪啪啪视频大全美女被插视频 国模掰阴露鲍艺术摄影 yinqi网 最新在线偷拍自拍在线观看 陌陌约炮视频在线观看 电影日本强奷免费播放 黯然销魂网站 爆乳操逼视频 saozixaoshuo sseyeyeshipin 超硪碰视频成人福利 QUFtYWduZXQ6P3h0PXVybjpidGloOjZFOUQ5NkVGODNCNzA4Nzc0RTExQzM1RkNFMEJDNTZC 步兵番号dky 日本销魂床上视频 爱人自拍种子 风流寡妇三级观看 欧美成人直播色图 腐男操逼视频 狠狠干迅雷下载 亚洲快播人与狗影视网 在线国产丝袜vr 在 线在日本一级黄片 白肤美性感情趣内衣美少妇和炮友直播啪啪啪 午夜伦理裸逼视频 国产调教手机在线 草地上的美脚视频在线播放 国产主播啪啪 爆乳中文字幕在线视频 无码 BT 乳交 种子 fulisisi mp4 爆乳波霸床上视频 五月天亚洲图片婷婷 幽冥鬼船激情片段 自拍 偷拍 亚洲 欧美 色 我已婚的妻子511视频 黄色片随意看高清版 动漫里番无修在线播放 东方成人α√ 粉鲍鱼d暗暗射动漫里番无修在线播放 东京热高清一本道免费 色逼爽片如何下载免费 萝莉白色丝袜h 美女被男人插屁视频app 萝莉成人在线视频 五月丁香蕉伊夜不色 漫漫自由路手机在线观看 美乳女神思瑞 黄色亚洲 八英里 影音先锋 藤井爱纱迅雷 西瓜影院福利导航 伦理片巨乳 360av免费 背德朋友的妈妈 伊色香蕉网 大桥未久生物教师 猫 亚洲色图 坏坏哥电影院 狠狠日 mp4 a熟女人妻av 丁香五月色拉 阅B影院 78ccdd 日本乱伦播放 熊黛林三级床视频戏 乐鱼 成人在线视频 FC2云播手机在线 色播久久爱视频精品 韩国少妇自慰在线电影 韩国主播福利视频禁片 19P农家炕上操进老娘屁眼 winktv荷恩视频在线 主播走光 磁力 大啪橹 威客理论片 福利电影不用下载插件 好屌操视颍 snis-573 在线观看 手机色视频日韩电影 福利king在线91视频 色大哥亚洲影院 石黑京香午夜在线视频 WWW_33ABAB_COM 21种女子美穴 插得好爽在线小说 佐佐木希最激情的作品 伸不进手塔用内视镜能看到各部位吗 9岁幼交实战 中年性感的妇女浪 制服丝袜电影 无码 肛交口交颜射人体艺术 亚洲 基地 婷婷 欧美性爱图10 第770弾み1 24p 我成了爸的性奴 我想看柳岩的屄 色肉体 女人秘密处动态照 骚情深爱色妹妹 WWW_168_INFO 老婆的三洞全被塞满 d81c19560000c525 淫淫网日日色黑人激情 爱与激情老二撸 迷奸恋足 明星 逼逼 俄罗斯嫰少女 美人妈妈和好色儿子转变为乱伦 狠狠鲁的小说婷五月 美女隐私图片欧美 老年男女qq聊性群号码 操逼网站来几个 欧美性爱15p写真 人妖变装图片 苍井空所有的三级片 俄罗斯美女 在线 京泉穴 徐静波 公司车抵税 大连航运职业技术学院 高清美女露逼图 sese666 幼幼开苞另类套图 人体艺术图片人体艺术图片大全 97xxuu 大槻响 欢媾 色小姐网盘 优酷操图和书 刚干的小嫩逼 操干妈自拍视频 60璐色图 大乳房大屁屁性爱影院 第四淫团 日本女同性爱性交 风间优美裸体影音先锋 巨乳3p之敏感新娘 高清掰开洞穴图片 五月激情五月婷婷 舔美女老师逼小说 WWW51BOBOCOM 林攻信作曲 色七七原图乱轮图片 成人小说网狠狠碰 美女芝艺术照 云播欧美a片 激情日色色色色 欧美亚洲乱伦迅雷 少妇被老外 哟女逼逼 抽插淫水初夜处女 成人大鸡少女篇吧 御村ありさオフィシャル美少女性交视频 快播幼苞少女 级片俄罗斯三 良家骚熟妇自拍 成人美女贴图区 谁有黄色网站bt下载 强奸幼幼的视频种子 做爱怎样的抽插女人喜欢 天天操逼网快播 徐莹阴道里塞的是什么东西 WWW2016ACOM 亚洲专区一本道 狼友做爱无码 苍井空三点影院 人猿泰山h版完整版影音先锋 147人体网 性吧有你自拍视频区 亚洲人体艺术诱惑 色色棒 激情淫乱色网站 插妹妹a片96偷拍美女上厕所 赤裸美臀坐脸 狼狗谢欣小说 内射人妻文字 成人网站幼女系列 2017丝袜少妇贴吧 日本女学生分享老师未经审查的阴茎 女人高潮男人狠狠射 熟女少妇淫乱综合片 撸管千百撸 av千百度 老外老头恋老图片 羽月馅丝空姐诱惑 亚洲最大的色的网wwwludashi147com 成人性av教育av网www8111rrcom 张莜雨性爱视频 抽插新婚女老师magnet 亚洲性va在线观看百度 姐妹骑 大波裸女视频 黄色片人体艺术摄影 丝袜诱惑狠狠射 巴欲影视天天影视 乳房性交乱伦小说 搜狗色妹妹 德国四级老片先锋 成人凌虐艳母小说 pa片资源吧 什么动作照相好看 姐姐色妹妹色骑姐骑妹 av妈妈乱 性动态图50p 搜索WWW0d3Pm1183020 抽插咪咪 sese尤物 寻抚摸女人乳房性爱视频 欧美美女劲爆人体 久久人居 爆操金发 夜夜直播抠b平台 美味的邂逅秒播影视 百度云盘片子分享 老公公操儿媳妇日本电影 黄色网站淫色狠狠 超大尺度裸体美女无水印 好大的奶好爽妇人 色中色影院最新资源免费观看 杭州阿里妈妈软件服务有限公司v 色乳爆帝国 欧美大鸡吧哥哥射射射 撸撸鸟性爱 AV网插妹妹www29azcom 日本人添逼逼视屏 av在线乐视影音 黑人强奸日本女人电影 www28secc 色导网肏妹妹 嫩逼骚妇的精致美逼特写 WWW_KKK_15_COM 功夫足球手机在线观看 花和尚com 5252she 欧美另类日韩 熟女撸撸色 AV丝在线 大鸡巴-百度 wee95sebacom 斗鱼三骚静静福利视频 涩影大全 wwwbu510com 手机影音先锋资源中文 1234zamagnet 色小组色天使狠狠干 少女小阴鲍p 妹妹给我乳交 女王调奴在线视频看 全本岳母系列肉小说 新cg儿 wwwhaole04com luo裸贷视频完整版 50p电脑版 wwwmaomiabcom 青青草免费视频在线观看好屌色 凤凰成人片 乱轮男女 公公操儿媳bb 美女爱操bxxaa552com 酥胸淫荡 欧美在线潮吹AV 日日性 青春草在线华人 3366se图片 青娱乐色在线 另类偷拍视频 木叶性爱图 欧洲美国色图 爱av陶宝 wwwzuixinAvcom 免费看顶级片 电影韩国天堂的卧室 精品套图奇米影视 淫色淫色永久久 淫深在线视频 av宗合 www111wo 亚洲天堂av无码免费 5XX 无码av动漫影音先锋 松岛枫AV手机观看 萝莉裸下体 泷泽萝拉第四部qvod 亚洲色图美腿丝袜欧洲无码 操小嫩穴视频 51cao超碰手机免费视频 色中色l 三级片赤裸的羔羊 亚洲色团30pwwwmcc676com 亚洲色图欧洲色图偷拍自拍另类色图在线 wwzzzzcfcf 日本女孩资源网站 汇聚全球 免费33qqrrcom日本 色艺无间正面是哪集 黄色激情漫画色女人网站激情小游戏 摩擦影院影音先锋色妹妹 肛交四姐妹下载 丰满中年骚熟妇尿尿 另类毛片网 巨乳孕妇av 444uh 520高清官网 淫荡人妻系列 互奸 kaixinjiqingxiaoshuo 大色妹电影网 大大小小色色色色色乱伦小说 爱操小视频 日b做爱吃鸡巴添bb 骚嫂子电影院 魔女栖息的家苍井空 偷拍自拍图片酒色 熟母检察官竹内さき快播电影 xfplay大长金 苍井空逼逼在线播放 亚洲色图欧美色吐激情小说 欧美人体艺术骚猛图 美女轮奸抽搐快播 淫淫亚洲色请网 幼交21先锋影音 youjizzmobilefreethunder 偷拍orrent 快播 逼逼逼逼逼逼逼 0809嫩逼导航 WWW_D4UU_COM 爱唯侦查成人小说 亚洲影院涩涩爱 撸叽av吉吉影音 爱舔黑木耳 爸爸让我屄妈妈 佐佐木遥中文电影 性交组图 yen道图片曾晓薇 狠狠捂逼 迷晕妹妹狂舔丝袜脚 日比四色鲁鲁网 操bi视频 乱伦娘 又能看的黄色网站部 少妞求操 xxx小次郎 日美女阴道 弟弟姐姐操逼小说 自拍偷拍 露脸 夫妻 www110bbcow 魅腿族 刂丶小说我和老师 大喰熟女食大便 什么样女人的屄最好看 东京热b图 幼幼幼女快播 老婆偷和外人肏屄 我爱肏小姨子的屄 狠撸撸蜜桃 母夜吟百度贴吧 鸡巴肏屄的图片 www韩国骚妇 老婆乳房不嫩 WWWSETUNET 美女被插内射p 第4色橹大妈 黑人鸡巴爆操 擂擂bt成人阁 摸了女儿的小屄 222bbbb最新网站是什么 裸模写真欧洲 乱伦成人网小说 成人图片插嫩穴 xiongmao29hd 李雨悠 五月天肉体 蛇妖人体艺人摄影 女人身体里面有狗精子 在线欧美性爱一级片 光棍电影影院yy俄罗斯年青母亲 彭守交中山 大帝av大帝在线视频成人 WWWNJHCZCCOM 3d肉葡萄团电影完整版电影影音先锋 幼木耳图片美女 妈妈爱鸡巴日 汤加丽完整版 wwavtt3orn 操妈妈屁眼小说 kancaobipian 国模莎莎大尺度 黑鸡巴曰大屄 18岁美女奶子图 大胆偷拍丰满少妇 韩国美女少妇人体艺术图片 亚州av无码mavtt155com 世界上最找的A片 美女制服诱惑男女 恋足thunderftp 性感的家庭私人女秘书 少女呻吟av domainwwwgulousecom 一线天av女优 2017偷拍自怕 久久在线视频精品99re6mweibocn 白嫩大腿熟妇照片 狂上师母 不用播放器的亚洲性爱 台湾黄色小说网 overthumbspenis图 结婚视频亚洲综合操的视频 鸡巴狂插无毛小穴 第二书包网站首页 亚洲成人18禁 制服丝袜挒表 色15p亚洲图片 wwwseqingwangzcom 傲娇淫荡少女 偷拍自拍美罗城 老电影连环画大全 欧美亚洲斯嘉丽 超碰日日嫂嫂 免费肉番伦理在线观看 av空姐制服电影 激情图片亚洲AV 一级强奸大片子 nuxingau 久久国产夫妻自拍刺激视频 乱伦片性感 美女少妇自慰舌舔自慰 罗伟春色 www.vipshop. com www淘宝网cn 色色在线 校园性爱 在线短片 六色女影院 美女人体网 四门电影网? 吉吉影音 伦理 亚洲 寂处阁 喜爱se在线播放 日本Av饭免费观看 上班不准穿内裤H文 超碰壁自拍 多多成人网 www001999C0M abp347图解 换妻自拍神马电影 偷窥自拍20 2828电影网金瓶艳史 99丫国产自拍 天天射天天日把b视频 下载 全网第一直播福利视频 小宫由梨? 在线福利导航 8ppav com 尼尔的福利视频的种子 rp社免费视频在线 汤姆影视三寸萝莉 双马尾萝莉揉逼视频 色老板免费线观看免费 天天啪啪,天天录,天天射,一本道。 藏经阁在线资源观看 日老妇女直播 日本医生鸭嘴器视频 紫薇阁精品自拍视频 猫咪无码资源 magnet 桃谷绘里香激情 日本AV黄图 mp4 9292看看午夜视频 h影院整片 影音先锋久久爱视频 ZUKO-100在线观看 婷婷六月 成人动漫 妖媚视频英文站 乱伦电影被女儿的男朋友抽插的妈妈 欧美性色黄大片 乌克兰a片 受美利坚共和国保护的av网站 亚洲强奸乱伦比比资源站 2017年色偷偷最新地址 天天影视av 国产自拍青青茶 免费视频人人橾123 国漫丝袜美腿 小仙儿西野翔 成人社区网站 66yxyx av10784 菲菲影院 二本道亚洲视频在线观看 电影风骚熟女诱惑 东方在线无码亚洲日韩 真实短片无码在线 日本人强奸美女的视频 小奶屄 gav成人小视频在线福利 2018仙女屋19禁电影大全 yy2525影院 第四色网男人香蕉 女捜査官波多野结衣mp4 色狗狗综合影视 午夜男人日b视频 色戒av电影 自拍摸BB在线观看 fuliwangzhan aiavatat av剃天堂 韩国美女主播金荷娜AV 韩国深夜综剧迅雷下载 成人短片库m4格式 福利自拍网 神马sm 超bt的日本摸奶节目 六间房系列 mp4 国产真实3p大战完整高清版本红衣 av视黄色视频 洋老外米糕新作品找不到小姐把酒店的服务员给草了 真人下体抽插下阴喷水视频 龚玥菲新金瓶百度云盘 乱片导行 日本zz╳x╳ 伦理站3166 www11peng av宝库 對來家庭寄宿的黑人巨屌發情了的母親 [中文字幕] 免费av站在线看2019未成年人禁止入内 小妹子福利视频 性感美女家中换衣服和土豪男友啪啪 性爽奴影院 性插怡红院视频 性感漂亮的美少妇驾照一直老板考不过与教练车震 小美女放学后被狂干 真实自拍做 正在播放 国模私拍 9|锅福利社 爆操视频 欧美肛交孕妇视频 讯雷av啪啪 做爱视频美女毛片 成人操逼视频app下载 日本中年a v 熟女 色戒网视频在线观看 色先锋 ckplayer 四虎图库 91顶级超模新作第242部第1集 av色琪琪20无码 被窝综合网站 操我求你了好痒用震动棒干我好爽飞了高潮了哦哦啊啊奥视频 插 女 自拍 国产 看米奇资源站 丁香wyt 哪部片子讲妻子被丈夫赶出家被邻居 thz-35中文 magnet www`脳脳X qiangjianyounvdianying wwse560info 操你啦日贺晓杨 槽榴社区 S*M、人*兽、奸*尸、幼*女、男*同 VENU无码 91超碰福利手机在线视频观看 迅雷XXOO视频种子 王瑞儿磁力链接 下载 女人自慰视频手机观看 久久操视频在线现观 拾红院 自慰78p 播放巨乳4 人狗性交电影网 男同事磨得我下面全是水视频 贞德h视频 午夜影院男人舔逼 超长阴茎的人妖视频3 一七六九b大香蕉在线 real-423 韩国美女主播杰西卡视频 用嘴就能让你爽 波多野结衣版在线播 宾馆情侣小视频福利 韩国美女激情vip秀187 无码 BT 乳交 种子 百度男人综合vA 自拍5678 YouasiaNi偷偷摸摸地去健身房WC 古装A片迅雷在线观看 青娱乐吧· 4438x8成人最新 五姑娘AV导航 人人色久久色 国产无码 在线 偷拍 捆绑天海翼 a片制服丝袜美腿av影院 女生7788被窝电影网_女生私人影院在线播放_66游戏网 西瓜影音 毛毛片国产自拍在线 美女白丝浴室诱惑视频 色吊丝国产在线视频 露b女主播磁力 美女被插福利 精品旅馆珍稀粉红房偷拍 精品网红在线 青青草人与兽人人干 日本无码福利自拍第一页 免费观看三级日本a片 国产3p视频ftp 日本efXXX Cosplay的h片 桃乃木香奈时间静止器 黄色视频软件速播影院 欧美手机伦理片 里番福利导航 国产av自拍在线 日本强奸制服丝袜电影 三级按摩推油伦理电影 黄色色播站 天降艳福不是福绝版在线看 堕落91社区 偷怕自拍激情小说 空中裸体视频 24p自在线拍 男女性交做爱插B吃鸡巴不堪的场面最污的 免费在线五码视频 国内十大禁片a片视频 国内中年夫妻自拍直接播放 国产自拍 主播诱惑’ 手机在线永久国产免费视频 成人四色快播 韩国直播 xxx 西瓜影音 fei插入式性交 747福利电影网 激情午夜夭 caoporn97齐鲁 豹纹内裤影院 性爱乱伦视频 欧美 av 无码 - Yahoo Search Results 三级片视频图片 女优在线论坛 6度影院鲁鲁片 人体写真美鲍 外国女人的比为什么总是没水 有声huangse小说 爱爱打炮影院 波多野结美图 韩国黄色强奸幼女 偷拍电影第一页 儿媳屄肉棒 白灵人体图片大胆 免费的淫秽图片 波多野结衣av无码百度云