It's _still_ not precise enough!

I spent a buttload of time working on the Tobler projection. It is now
unbearably slow and still has an average distortion of, like .04. It's
okay, though. I've got a plan to fix both of those problems. It involves
ordinary differential equations and inheritance. Nyeh heh heh heh heh!

src/apps/MapAnalyzer.java

@@ -75,11 +75,11 @@ public class MapAnalyzer extends MapApplication {
 	private static final int CHART_WIDTH = 400;
 	private static final int ROUGH_SAMP_NUM = 250;
 	private static final int FINE_SAMP_NUM = 1000;
-	private static final double GLOBE_RES = .02;
+	private static final double GLOBE_RES = .01;
 	
 	private static final FileChooser.ExtensionFilter[] RASTER_TYPES = {
-			new FileChooser.ExtensionFilter("PNG", "*.png"),
+			new FileChooser.ExtensionFilter("JPG", "*.jpg","*.jpeg","*.jpe","*.jfif"),
-			new FileChooser.ExtensionFilter("JPG", "*.jpg","*.jpeg","*.jpe","*.jfif") };
+			new FileChooser.ExtensionFilter("PNG", "*.png") };
 	
 	private static final Projection[] PROJ_ARR = { Projection.MERCATOR, Projection.PLATE_CARREE, Projection.GALL_PETERS,
 			Projection.HOBO_DYER, Projection.BEHRMANN, Projection.LAMBERT_CYLIND, Projection.E_A_CYLIND,
@@ -122,9 +122,9 @@ public class MapAnalyzer extends MapApplication {
 		this.updateBtn = buildUpdateButton(this::calculateAndUpdate);
 		this.updateBtn.setText("Calculate"); //I don't need to follow your darn conventions!
 		final Button saveMapBtn = buildSaveButton(true, "map", RASTER_TYPES,
-				Procedure.NONE, this::calculateAndSaveMap);
+				RASTER_TYPES[1], Procedure.NONE, this::calculateAndSaveMap);
 		final Button savePltBtn = buildSaveButton(true, "plots", RASTER_TYPES,
-				Procedure.NONE, this::calculateAndSavePlot);
+				RASTER_TYPES[1], Procedure.NONE, this::calculateAndSavePlot);
 		final HBox buttons = new HBox(5, updateBtn, saveMapBtn, savePltBtn);
 		buttons.setAlignment(Pos.CENTER);
 		
@@ -212,7 +212,7 @@ public class MapAnalyzer extends MapApplication {
 			sizeChart.getData().add(histogram(distortionG[0],
 					-2, 2, 14, Math::exp));
 			shapeChart.getData().add(histogram(distortionG[1],
-					0, 2.8, 14, (x) -> 1/(x*x/2+1-x*Math.sqrt(x*x/4+1))));
+					0, 2, 14, (x) -> 1/(x*x/2+1-x*Math.sqrt(x*x/4+1))));
 			
 			avgSizeDistort.setText(format(Math2.stdDev(distortionG[0])));
 			avgShapeDistort.setText(format(Math2.mean(distortionG[1])));
@@ -275,14 +275,14 @@ public class MapAnalyzer extends MapApplication {
 				
 				final int r, g, b;
 				if (sizeDistort < 0) { //if compressing
-					r = (int)(255.9*Math.exp(-shapeDistort/1.5));
+					r = (int)(255.9*Math.exp(-shapeDistort*.6));
-					g = (int)(255.9*Math.exp(-shapeDistort/1.5)*Math.exp(sizeDistort/1.5));
+					g = (int)(255.9*Math.exp(-shapeDistort*.6)*Math.exp(sizeDistort*.6));
 					b = g;
 				}
 				else { //if dilating
-					r = (int)(255.9*Math.exp(-shapeDistort/1.5)*Math.exp(-sizeDistort/1.5));
+					r = (int)(255.9*Math.exp(-shapeDistort*.6)*Math.exp(-sizeDistort*.6));
-					g = r;
+					g = r; //I find .6 to ba a rather visually pleasing sensitivity
-					b = (int)(255.9*Math.exp(-shapeDistort/1.5));
+					b = (int)(255.9*Math.exp(-shapeDistort*.6));
 				}
 				
 				final int argb = ((((((0xFF)<<8)+r)<<8)+g)<<8)+b;

src/apps/MapApplication.java

@@ -64,8 +64,8 @@ public abstract class MapApplication extends Application {
 	private static final String[] ASPECT_NAMES = { "Standard", "Transverse", "Center of Mass", "Jerusalem", "Point Nemo",
 			"Longest Line", "Longest Line Transverse", "Cylindrical", "Conic", "Tetrahedral", "Quincuncial", "Antipode", "Random" };
 	private static final double[][] ASPECT_VALS = {
-			{ 90., 0., 29.9792, 31.7833, 48.8767, -28.5217,-46.4883,-35.,    -10.,  47., 60. },
+			{ 90., 0., 29.98, 31.78, 48.88, -28.52,-46.4883,-35.,  -10.,  47., 60. },
-			{  0., 0., 31.1344, 35.2160, 56.6067, 141.451,  16.5305,-13.6064, 65.,-173., -6. },
+			{  0., 0., 31.13, 35.22, 56.61, 141.45, 16.5305,-13.61, 65.,-173., -6. },
 			{  0., 0.,-32.,  -35.,  -45.,   161.5, 137.,    145., -150., 138.,-10. } };
 	
 	
@@ -107,6 +107,7 @@ public abstract class MapApplication extends Application {
 	 */
 	protected Region buildInputSelector(
 			FileChooser.ExtensionFilter[] allowedExtensions,
+			FileChooser.ExtensionFilter defaultExtension,
 			Consumer<File> inputSetter) {
 		final Label label = new Label("Current input:");
 		final Text inputLabel = new Text("None");
@@ -115,6 +116,7 @@ public abstract class MapApplication extends Application {
 		inputChooser.setInitialDirectory(new File("input"));
 		inputChooser.setTitle("Choose an input map");
 		inputChooser.getExtensionFilters().addAll(allowedExtensions);
+		inputChooser.setSelectedExtensionFilter(defaultExtension);
 		
 		final Button loadButton = new Button("Choose input...");
 		loadButton.setTooltip(new Tooltip(
@@ -300,13 +302,15 @@ public abstract class MapApplication extends Application {
 	 */
 	protected Button buildSaveButton(boolean bindCtrlS, String savee,
 			FileChooser.ExtensionFilter[] allowedExtensions,
+			FileChooser.ExtensionFilter defaultExtension,
 			Procedure settingCollector,
 			BiConsumer<File, ProgressBarDialog> calculateAndSaver) {
 		final FileChooser saver = new FileChooser();
 		saver.setInitialDirectory(new File("output"));
-		saver.setInitialFileName("my"+savee+allowedExtensions[0].getExtensions().get(0).substring(1));
+		saver.setInitialFileName("my"+savee+defaultExtension.getExtensions().get(0).substring(1));
 		saver.setTitle("Save "+savee);
 		saver.getExtensionFilters().addAll(allowedExtensions);
+		saver.setSelectedExtensionFilter(defaultExtension);
 		
 		final Button saveButton = new Button("Save "+savee+"...");
 		saveButton.setOnAction(new EventHandler<ActionEvent>() {

src/apps/MapDesignerRaster.java

@@ -63,8 +63,8 @@ public class MapDesignerRaster extends MapApplication {
 	
 	
 	private static final FileChooser.ExtensionFilter[] RASTER_TYPES = {
-			new FileChooser.ExtensionFilter("PNG", "*.png"),
+			new FileChooser.ExtensionFilter("JPG", "*.jpg","*.jpeg","*.jpe","*.jfif"),
-			new FileChooser.ExtensionFilter("JPG", "*.jpg","*.jpeg","*.jpe","*.jfif") };
+			new FileChooser.ExtensionFilter("PNG", "*.png") };
 	
 	private static final Projection[] PROJ_ARR = { Projection.MERCATOR, Projection.PLATE_CARREE, Projection.HOBO_DYER,
 			Projection.GALL, Projection.STEREOGRAPHIC, Projection.POLAR, Projection.E_A_AZIMUTH,
@@ -102,13 +102,16 @@ public class MapDesignerRaster extends MapApplication {
 	@Override
 	protected Node makeWidgets() {
 		this.aspect = new double[3];
-		final Node inputSelector = buildInputSelector(RASTER_TYPES, this::setInput);
+		final Node inputSelector = buildInputSelector(RASTER_TYPES,
-		final Node projectionSelector = buildProjectionSelector(PROJ_ARR, Projection.MERCATOR, Procedure.NONE);
+				RASTER_TYPES[0], this::setInput);
-		final Node aspectSelector = buildAspectSelector(this.aspect, Procedure.NONE);
+		final Node projectionSelector = buildProjectionSelector(PROJ_ARR,
+				Projection.MERCATOR, Procedure.NONE);
+		final Node aspectSelector = buildAspectSelector(this.aspect,
+				Procedure.NONE);
 		final Node parameterSelector = buildParameterSelector(Procedure.NONE);
 		this.updateBtn = buildUpdateButton(this::updateMap);
 		this.saveMapBtn = buildSaveButton(true, "map", RASTER_TYPES,
-				this::collectFinalSettings, this::calculateAndSaveMap);
+				RASTER_TYPES[1], this::collectFinalSettings, this::calculateAndSaveMap);
 		final HBox buttons = new HBox(5, updateBtn, saveMapBtn);
 		buttons.setAlignment(Pos.CENTER);
 		
@@ -151,7 +154,8 @@ public class MapDesignerRaster extends MapApplication {
 	
 	
 	private void updateMap() {
-		display.setImage(map());
+		display.setImage(makeImage(
+				this.getProjection().map(IMG_WIDTH, this.getParams(), aspect)));
 	}
 	
 	
@@ -163,9 +167,13 @@ public class MapDesignerRaster extends MapApplication {
 	
 	
 	protected void calculateAndSaveMap(File file, ProgressBarDialog pBar) {
-		Image theMap = map(
+		final int[] outDims = configDialog.getDims();
-				configDialog.getDims(), configDialog.getSmoothing(), pBar); //calculate
+		final int smoothing = configDialog.getSmoothing();
+		double[][][] points =this.getProjection().map(
+				outDims[0]*smoothing, outDims[1]*smoothing, this.getParams(), aspect, pBar::setProgress);
 		pBar.setProgress(-1);
+		Image theMap = makeImage(points, smoothing); //calculate
+		
 		final String filename = file.getName();
 		final String extension = filename.substring(filename.lastIndexOf('.')+1);
 		try {
@@ -179,66 +187,42 @@ public class MapDesignerRaster extends MapApplication {
 	}
 	
 	
-	public Image map() {
+	private Image makeImage(double[][][] points) {
-		final double a = this.getProjection().getAspectRatio(this.getParams());
+		return makeImage(points, 1);
-		return map(new int[] { IMG_WIDTH, (int)(IMG_WIDTH/a) }, 1);
 	}
 	
-	public Image map(int[] outputDims, int smoothing) {
+	private Image makeImage(double[][][] points, int step) {
-		return map(outputDims,smoothing, null);
+		final PixelReader in = input.getPixelReader();
-	}
+		final WritableImage out = new WritableImage(points[0].length/step, points.length/step);
-	
+		for (int y = 0; y < out.getHeight(); y ++) {
-	public Image map(int[] outDims, int smoothing,
+			for (int x = 0; x < out.getWidth(); x ++) {
-			ProgressBarDialog pbar) {
+				int[] colors = new int[step*step];
-		final Projection proj = this.getProjection();
+				for (int dy = 0; dy < step; dy ++) {
-		final double[] params = this.getParams();
+					for (int dx = 0; dx < step; dx ++) {
-		final PixelReader ref = input.getPixelReader();
+						colors[step*dy+dx] = getArgb(points[step*y+dy][step*x+dx], in, input.getWidth(), input.getHeight());
-		final int[] refDims = {(int)input.getWidth(), (int)input.getHeight()};
-		
-		WritableImage img = new WritableImage(outDims[0], outDims[1]);
-		
-		for (int x = 0; x < outDims[0]; x ++) {
-			for (int y = 0; y < outDims[1]; y ++) {
-				int[] colors = new int[smoothing*smoothing];
-				int i = 0;
-				for (double dx = 0.5/smoothing; dx < 1; dx += 1.0/smoothing) {
-					for (double dy = .5/smoothing; dy < 1; dy += 1.0/smoothing) {
-						colors[i] = getArgb(x+dx, y+dy,
-								proj,params,aspect,ref,refDims,outDims);
-						i ++;
 					}
 				}
-				img.getPixelWriter().setArgb(x, y, blend(colors));
+				out.getPixelWriter().setArgb(x, y, blend(colors));
 			}
-			if (pbar != null)	pbar.setProgress((double)(x+1)/outDims[0]);
 		}
-		
+		return out;
-		return img;
 	}
 	
 	
-	public static int getArgb(double x, double y, Projection proj, double[] params, double[] pole,
+	public static int getArgb(double[] coords, PixelReader in,
-			PixelReader ref, int[] refDims, int[] outDims) {
+			double inWidth, double inHeight) { // returns the color of any coordinate on earth
-		final double[] coords = proj.inverse(
+		if (coords == null) 	return 0;
-				2.*x/outDims[0]-1, 1-2.*y/outDims[1], params, pole);
-		if (coords != null)
-			return getColorAt(coords, ref, refDims);
-		else
-			return 0;
-	}
 		
-	
-	public static int getColorAt(double[] coords, PixelReader ref, int[] refDims) { // returns the color of any coordinate on earth
 		double x = 1/2.0 + coords[1]/(2*Math.PI);
-		x = (x - Math.floor(x)) * refDims[0];
+		x = (x - Math.floor(x)) * inWidth;
 		
-		double y = refDims[1]/2.0 - coords[0]*refDims[1]/(Math.PI);
+		double y = inHeight/2.0 - coords[0]*inHeight/(Math.PI);
 		if (y < 0)
 			y = 0;
-		else if (y >= refDims[1])
+		else if (y >= inHeight)
-			y = refDims[1] - 1;
+			y = inHeight - 1;
 		
-		return ref.getArgb((int)x, (int)y);
+		return in.getArgb((int) x, (int) y);
 	}
 	
 	

src/apps/MapDesignerVector.java

@@ -104,12 +104,15 @@ public class MapDesignerVector extends MapApplication {
 	@Override
 	public Node makeWidgets() {
 		this.aspect = new double[3];
-		final Node inputSelector = buildInputSelector(VECTOR_TYPES, this::setInput);
+		final Node inputSelector = buildInputSelector(VECTOR_TYPES,
-		final Node projectionSelector = buildProjectionSelector(PROJ_ARR, Projection.MERCATOR, this::updateMap);
+				VECTOR_TYPES[0], this::setInput);
-		final Node aspectSelector = buildAspectSelector(this.aspect, this::updateMap);
+		final Node projectionSelector = buildProjectionSelector(PROJ_ARR,
+				Projection.MERCATOR, this::updateMap);
+		final Node aspectSelector = buildAspectSelector(this.aspect,
+				this::updateMap);
 		final Node parameterSelector = buildParameterSelector(this::updateMap);
 		this.saveBtn = buildSaveButton(true, "map", VECTOR_TYPES,
-				Procedure.NONE, this::calculateAndSaveMap);
+				VECTOR_TYPES[0], Procedure.NONE, this::calculateAndSaveMap);
 		
 		final VBox layout = new VBox(5,
 				inputSelector, new Separator(), projectionSelector,
@@ -220,12 +223,18 @@ public class MapDesignerVector extends MapApplication {
 		g.clearRect(0, 0, c.getWidth(), c.getHeight());
 		for (List<double[]> closedCurve: img) {
 			g.beginPath();
-			for (int i = 0; i < closedCurve.size(); i ++) {
+			g.moveTo(closedCurve.get(0)[0], closedCurve.get(0)[1]);
-				double[] p = closedCurve.get(i);
+			for (int i = 1; i < closedCurve.size()+1; i ++) {
-				if (i == 0)	g.moveTo(p[0], p[1]);
+				double[] p0 = closedCurve.get(i-1);
-				else		g.lineTo(p[0], p[1]);
+				double[] p1 = closedCurve.get(i%closedCurve.size());
+				if (Math.hypot(p1[0]-p0[0], p1[1]-p0[1]) >= c.getWidth()/10) {
+					g.stroke();
+					g.beginPath();
+					g.moveTo(p1[0], p1[1]);
+				}
+				else
+					g.lineTo(p1[0], p1[1]);
 			}
-			g.closePath();
 			g.stroke();
 		}
 	}

src/maps/Projection.java

@@ -25,6 +25,7 @@ package maps;
 
 import java.util.ArrayList;
 import java.util.List;
+import java.util.function.DoubleConsumer;
 import java.util.function.UnaryOperator;
 
 import org.apache.commons.math3.complex.Complex;
@@ -109,8 +110,8 @@ public enum Projection {
 	},
 	
 	E_A_CYLIND("Equal-area cylindrical", "A generalized equal-area cylindrical projection",
-			Math.PI, 0b1111, "cylindrical", "equal-area",
+			1.977, 0b1111, "cylindrical", "equal-area",
-			new String[]{"Std. parallel"}, new double[][]{{0, 75, 37.5}}) {
+			new String[]{"Std. parallel"}, new double[][]{{0, 85, 37.5}}) {
 		public double[] project(double lat, double lon, double[] params) {
 			final double a = Math.pow(Math.cos(Math.toRadians(params[0])), 2);
 			return new double[] {lon, Math.sin(lat)/a};
@@ -395,17 +396,41 @@ public enum Projection {
 	},
 	
 	TOBLER("Tobler", "An equal-area projection shaped like a hyperellipse",
-			2., 0b1001, "pseudocylindrical", "equal-area",new String[]{"alpha","gamma","k"},
+			2., 0b1001, "pseudocylindrical", "equal-area",new String[]{"alpha","K"},
-			new double[][] {{0,1,0}, {1,5,2.5}, {1,2,1.831}}) {
+			new double[][] {{0,1,0}, {1,5,2.5}}) {
 				public double[] project(double lat, double lon, double[] params) {
-					return new double[] {
+					final double g = 1/NumericalAnalysis.simpsonIntegrate(0, 1,
-							lon*Tobler.xfacFromLat(lat)*18,
+							Tobler::hyperEllipse, .01, params);//I think my gammaay be defined differently from Tobler's
-							Tobler.yFromLat(lat)*Math.PI/10 };
+					final double y =
+							NumericalAnalysis.newtonRaphsonApproximation(
+									Math.abs(Math.sin(lat)), Math.abs(Math.sin(lat)),
+									Tobler::sinPhi, Tobler::dsinPhidY, .01,
+									params[0], params[1], g)*Math.signum(lat);
+					return new double[] { Tobler.X(y, lon, params), Math.PI/2*y }; //I suppose I could make this about twice as fast if I really wanted to deal with more inheritance, but this is fast enough.
 				}
 				public double[] inverse(double x, double y, double[] params) {
-					return new double[] {
+					return new double[] { 0, Tobler.lam(x,y,params) };
-							Tobler.latFromY(5*y), //TODO: make this be real
+				}
-							x/Tobler.xfacFromY(5*y)*Math.PI/18 };
+				public double[][][] map(double w, double h, double[] params, double[] pole, DoubleConsumer tracker) { //generate a matrix of coordinates based on a map projection
+					final int n = (int) h; //just so I don't forget
+					final double g = 1/NumericalAnalysis.simpsonIntegrate(0, 1,
+							Tobler::hyperEllipse, .01, params);
+					double[] z = NumericalAnalysis.simpsonODESolve(
+							1, n, Tobler::dZdY, Math.min(1./n,.01),
+							params[0], params[1], g);
+					double[][][] output = new double[(int) h][(int) w][2]; //the coordinate matrix
+					for (int y = 0; y < output.length; y ++) {
+						final double zoy;
+						if (y < output.length/2) 	zoy = z[n - 2*y - 1];
+						else 						zoy =-z[2*y + 1 - n];
+						for (int x = 0; x < output[y].length; x ++) {
+							output[y][x] = inverse((2*x+1)/w-1, 1-(2*y+1)/h, params, pole);
+							output[y][x][0] = Math.asin(zoy);
+						}
+						if (tracker != null) 	tracker.accept((double) y/output.length);
+					}
+					
+					return output;
 				}
 	},
 	
@@ -469,9 +494,10 @@ public enum Projection {
 			final double cosphi0 = Math.cos(Math.toRadians(params[0]));
 			return NumericalAnalysis.newtonRaphsonApproximation(
 					x*Math.PI*(1 + cosphi0), y*Math.PI,
-					WinkelTripel::f1pX, WinkelTripel::f2pY, WinkelTripel::df1dphi,
+					y*Math.PI/2, x*Math.PI*(1 + Math.cos(y*Math.PI/2))/2,
-					WinkelTripel::df1dlam, WinkelTripel::df2dphi, WinkelTripel::df2dlam,
+					WinkelTripel::f1pX, WinkelTripel::f2pY,
-					.05, cosphi0);
+					WinkelTripel::df1dphi, WinkelTripel::df1dlam,
+					WinkelTripel::df2dphi, WinkelTripel::df2dlam, .05, cosphi0);
 		}
 		public double getAspectRatio(double[] params) {
 			return 1 + Math.cos(Math.toRadians(params[0]));
@@ -766,13 +792,26 @@ public enum Projection {
 	}
 	
 	
-	public double[][][] map(int size, double[] params) { //generate a matrix of coordinates based on a map projection
+	public double[][][] map(int size, double[] params) {
-		final int w = size, h = (int)(size/this.getAspectRatio(params));
+		return map(size, params, new double[] {Math.PI/2,0,0});
-		double[][][] output = new double[h][w][2]; //the coordinate matrix
+	}
 	
-		for (int y = 0; y < output.length; y ++)
+	public double[][][] map(int size, double[] params, double[] pole) {
+		final double ratio = this.getAspectRatio(params);
+		if (ratio < 1)
+			return map(Math.round(size*ratio), size, params, pole, null);
+		else
+			return map(size, Math.round(size/ratio), params, pole, null);
+	}
+	
+	public double[][][] map(double w, double h, double[] params, double[] pole, DoubleConsumer tracker) { //generate a matrix of coordinates based on a map projection
+		double[][][] output = new double[(int) h][(int) w][2]; //the coordinate matrix
+		
+		for (int y = 0; y < output.length; y ++) {
 			for (int x = 0; x < output[y].length; x ++)
-				output[y][x] = inverse(2*(x+.5)/w - 1, 1 - 2*(y+.5)/h, params); //s0 is this point on the sphere
+				output[y][x] = inverse((2*x+1)/w-1, 1-(2*y+1)/h, params, pole);
+			if (tracker != null) 	tracker.accept((double) y/output.length);
+		}
 		
 		return output;
 	}

src/maps/Tobler.java

@@ -23,7 +23,7 @@
  */
 package maps;
 
-import org.apache.commons.math3.analysis.interpolation.SplineInterpolator;
+import util.NumericalAnalysis;
 
 /**
  * A class of values and functions used to approximate the Tobler projection
@@ -32,36 +32,43 @@ import org.apache.commons.math3.analysis.interpolation.SplineInterpolator;
  */
 public class Tobler {
 
-	private static final double[][] TABLE = {
+	public static final double lam(double x, double y, double[] params) {
-			{  -90,-85,-80,-75,-70,-65,-60,-55,-50,-45,-40,-35,-30,-25,-20,-15,-10,-05, 00,
+		final double a = params[0];
-					05, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90 },
+		return Math.PI * x / (a + (1-a)*hyperEllipse(y, params));
-			{-5.000000,-4.921966,-4.782840,-4.610387,-4.410173,-4.184957,-3.940243,-3.678861,-3.398039,-3.103326,-2.794726,-2.475016,-2.144195,-1.802264,-1.454780,-1.096182,-0.734807,-0.367979, 0.000000,
-						0.367979, 0.734807, 1.096182, 1.454780, 1.802264, 2.144195, 2.475016, 2.794726, 3.103326, 3.398039, 3.678861, 3.940243, 4.184957, 4.410173, 4.610387, 4.782840, 4.921966, 5.000000 },
-			{ 0.000000, 0.022944, 0.030411, 0.035519, 0.039413, 0.042524, 0.045070, 0.047181, 0.048940, 0.050406, 0.051625, 0.052627, 0.053442, 0.054091, 0.054595, 0.054973, 0.055242, 0.055429, 0.055555,
-						0.055429, 0.055242, 0.054973, 0.054595, 0.054091, 0.053442, 0.052627, 0.051625, 0.050406, 0.048940, 0.047181, 0.045070, 0.042524, 0.039413, 0.035519, 0.030411, 0.022944, 0.000000 }
-		}; //these values were calculated with some kind of iterative approach in 1973.
-	
-	private static final SplineInterpolator SI = new SplineInterpolator(); //I don't understand why these can't be static functions, but fine
-	
-	
-	
-	public static final double xfacFromLat(double lat) {
-		return SI.interpolate(TABLE[0], TABLE[2]).value(Math.toDegrees(lat)); //I also don't understand why I can't call interpolate into a final Object instead of calling it every time. It's weird. If I try to do that, then invoking this class terminates the thread, like there's something so awful about calling that method outside of any method that Java just dies
 	}
 	
 	
-	public static final double yFromLat(double lat) {
+	public static final double X(double y, double lam, double[] params) {
-		return SI.interpolate(TABLE[0], TABLE[1]).value(Math.toDegrees(lat));
+		final double a = params[0];
+		return lam * (a + (1-a)*hyperEllipse(y, params));
 	}
 	
 	
-	public static final double latFromY(double pdfe) {
+	public static final double sinPhi(double y, double[] params) { //subtract sin(phi), once you have phi, to get error
-		return Math.toRadians(SI.interpolate(TABLE[1], TABLE[0]).value(pdfe));
+		final double a = params[0], g = params[2];
+		return (a*y + (1-a)*NumericalAnalysis.simpsonIntegrate(
+				0, y, Tobler::hyperEllipse, .0625, params))/
+				(a + (1-a)/g);
 	}
 	
 	
-	public static final double xfacFromY(double pdfe) {
+	public static final double dsinPhidY(double y, double[] params) {
-		return SI.interpolate(TABLE[1], TABLE[2]).value(pdfe);
+		final double a = params[0], g = params[2];
+		return (a + (1-a)*hyperEllipse(y, params))/
+				(a + (1-a)/g);
+	}
+	
+	
+	public static final double dZdY(double y, double[] params) {
+		final double a = params[0], g = params[2];
+		return (a + (1-a)*hyperEllipse(y, params))/
+				(a + (1-a)/g);
+	}
+	
+	
+	public static final double hyperEllipse(double y, double[] params) {
+		final double k = params[1];
+		return Math.pow(1 - Math.pow(Math.abs(y),k), 1/k);
 	}
 
 }

src/util/NumericalAnalysis.java

@@ -32,43 +32,123 @@ import java.util.Arrays;
  */
 public class NumericalAnalysis {
 
+	public static final void main(String[] args) {
+		System.out.println(simpsonIntegrate(-1,1, (x,prms) -> 4*Math.sqrt(1-x*x), .001, null));
+	}
 	/**
-	 * Applies Newton's method in two dimensions to solve for phi and lam given
+	 * Performs a definite integral using Simpson's rule and a constant step size
-	 * desired x and y values, x(phi,lam), y(phi,lam), and some derivatives
+	 * @param a The start of the integration region
-	 * @param x0 The x value that the functions need to match
+	 * @param b The end of the integration region (must be greater than a)
-	 * @param y0 The y value that the functions need to match
+	 * @param f The integrand
-	 * @param f1pX x in terms of phi and lam
+	 * @param h The step size (must be positive)
-	 * @param f2pY y in terms of phi and lam
+	 * @param constants Constant parameters for the function
+	 * @return \int_a^b \! f(x) \, \mathrm{d}x
+	 */
+	public static final double simpsonIntegrate(double a, double b, ScalarFunction f, double h, double... constants) {
+		double sum = 0;
+		for (double x = a; x < b; x += h) {
+			if (x+h > b) 	h = b-x;
+			sum += h/6*(f.evaluate(x,constants)
+					+ 4*f.evaluate(x+h/2, constants)
+					+   f.evaluate(x+h, constants));
+		}
+		return sum;
+	}
+	
+	
+	/**
+	 * Solves a simple ODE using Simpson's rule and a constant step size
+	 * @param T The maximum time value at which to sample (must be positive)
+	 * @param n The desired number of spaces (or the number of samples minus 1)
+	 * @param f The derivative of y with respect to time
+	 * @param h The internal step size (must be positive)
+	 * @param constants Constant parameters for the function
+	 * @return the double[] y, where y[i] is the value of y at t=i*T/n
+	 */
+	public static final double[] simpsonODESolve(double T, int n, ScalarFunction f, double h, double... constants) {
+		final double[] y = new double[n+1]; //the output
+		double t = 0;
+		double sum = 0;
+		for (int i = 0; i <= n; i ++) {
+			while (t < i*T/n) {
+				final double tph = Math.min(t+h, i*T/n);
+				sum += (tph-t)/6*(f.evaluate(t, constants)
+							  + 4*f.evaluate((t+tph)/2, constants)
+							  +   f.evaluate(tph, constants));
+				t = tph;
+			}
+			y[i] = sum;
+		}
+		return y;
+	}
+	
+	
+	/**
+	 * Applies Newton's method in one dimension to solve for x such that f(x)=y
+	 * @param y Desired value for f
+	 * @param x0 Initial guess for x
+	 * @param f The error in terms of x
+	 * @param dfdx The derivative of f with respect to x
+	 * @param tolerance The maximum error that this can return
+	 * @param constants Constant parameters for the function
+	 * @return The value of x that puts f near 0, or NaN if it does not
+	 * 		converge in 5 iterations
+	 */
+	public static final double newtonRaphsonApproximation(
+			double y, double x0, ScalarFunction f, ScalarFunction dfdx,
+			double tolerance, double... constants) {
+		double x = x0;
+		double error = Math.PI;
+		for (int i = 0; i < 5 && error > tolerance; i ++) {
+			error = f.evaluate(x, constants) - y;
+			final double dydx = dfdx.evaluate(x, constants);
+			x -= error/dydx;
+		}
+		if (error > tolerance)
+			return Double.NaN;
+		else
+			return x;
+	}
+	
+	
+	/**
+	 * Applies Newton's method in two dimensions to solve for phi and lam such
+	 * that f1(phi,lam)=x and f2(phi,lam)=y
+	 * @param x Desired value for f1
+	 * @param y Desired value for f2
+	 * @param phi0 Initial guess for phi
+	 * @param lam0 Initial guess for lam
+	 * @param f1 x-error in terms of phi and lam
+	 * @param f2 y-error in terms of phi and lam
 	 * @param df1dp The partial derivative of x with respect to phi
 	 * @param df1dl The partial derivative of x with respect to lam
 	 * @param df2dp The partial derivative of y with respect to phi
 	 * @param df2dl The partial derivative of y with respect to lam
 	 * @param tolerance The maximum error that this can return
-	 * @param params Constant parameters for the functions
+	 * @param constants Constant parameters for the functions
-	 * @return the values of phi and lam that put f1pX and f2pY near x0 and y0
+	 * @return the values of phi and lam that put f1 and f2 near 0, or
+	 * 			<code>null</code> if it does not converge in 5 iterations.
 	 */
-	public static final double[] newtonRaphsonApproximation(double x0, double y0,
+	public static final double[] newtonRaphsonApproximation(double x, double y,
-			VectorFunction fxpX, VectorFunction fypY, VectorFunction dfxdp,
+			double phi0, double lam0, VectorFunction f1, VectorFunction f2,
-			VectorFunction dfxdl, VectorFunction dfydp, VectorFunction dfydl,
+			VectorFunction df1dp, VectorFunction df1dl, VectorFunction df2dp,
-			double tolerance, double... params) {
+			VectorFunction df2dl, double tolerance, double... constants) {
-		double x = x0;
+		double phi = phi0;
-		double y = y0;
+		double lam = lam0;
-		double phi = y/2;
-		double lam = x/2; // I used equirectangular for my initial guess
 		double error = Math.PI;
 		
 		for (int i = 0; i < 5 && error > tolerance; i++) {
-			final double f1 = fxpX.evaluate(phi, lam, params) - x;
+			final double F1mx = f1.evaluate(phi, lam, constants) - x;
-			final double f2 = fypY.evaluate(phi, lam, params) - y;
+			final double F2my = f2.evaluate(phi, lam, constants) - y;
-			final double df1dP = dfxdp.evaluate(phi, lam, params);
+			final double dF1dP = df1dp.evaluate(phi, lam, constants);
-			final double df1dL = dfxdl.evaluate(phi, lam, params);
+			final double dF1dL = df1dl.evaluate(phi, lam, constants);
-			final double df2dP = dfydp.evaluate(phi, lam, params);
+			final double dF2dP = df2dp.evaluate(phi, lam, constants);
-			final double df2dL = dfydl.evaluate(phi, lam, params);
+			final double dF2dL = df2dl.evaluate(phi, lam, constants);
 			
-			phi -= (f1*df2dL - f2*df1dL) / (df1dP*df2dL - df2dP*df1dL);
+			phi -= (F1mx*dF2dL - F2my*dF1dL) / (dF1dP*dF2dL - dF2dP*dF1dL);
-			lam -= (f2*df1dP - f1*df2dP) / (df1dP*df2dL - df2dP*df1dL);
+			lam -= (F2my*dF1dP - F1mx*dF2dP) / (dF1dP*dF2dL - dF2dP*dF1dL);
 			
-			error = Math.hypot(f1, f2);
+			error = Math.hypot(F1mx, F2my);
 		}
 		
 		if (error > tolerance) // if it aborted due to timeout
@@ -94,8 +174,8 @@ public class NumericalAnalysis {
 	 * @param x The input value
 	 * @param X The sorted array of inputs on which to interpolate
 	 * @param f The sorted array of outputs on which to interpolate
-	 * @param from The index of the arrays at which to start
+	 * @param from The index of the arrays at which to start (inclusive)
-	 * @param to The index of the arrays at which to stop
+	 * @param to The index of the arrays at which to stop (exclusive)
 	 * @return f(x), approximately
 	 */
 	public static final double aitkenInterpolate(double x,
@@ -119,6 +199,11 @@ public class NumericalAnalysis {
 	
 	
 	
+	@FunctionalInterface
+	public interface ScalarFunction {
+		public double evaluate(double x, double[] constants);
+	}
+	
 	@FunctionalInterface
 	public interface VectorFunction {
 		public double evaluate(double x, double y, double[] constants);