Huygens-Fresnel wavefront tracing to calculate diffraction effects in telescopes.

CcdToImage.java

package com.clarke.agnes.experiments;

import javax.imageio.ImageIO;
import java.awt.image.BufferedImage;
import java.awt.image.WritableRaster;
import java.io.File;
import java.io.IOException;

public class CcdToImage {

    private int counter = 1;

    public void noiseToImage(ImageFile[] files, TelescopeSimulator.Photosite[][] ccd, double max, double min, String formatName, File outputFolder) {
        try {
            int dimension1 = ccd.length;
            int dimension2 = ccd[0].length;
            for (ImageFile file : files) {
                BufferedImage bi = new BufferedImage(dimension1, dimension2, BufferedImage.TYPE_3BYTE_BGR);
                WritableRaster wr = bi.getRaster().createCompatibleWritableRaster();
                file.setBufferedImage(bi);
                file.setWritableRaster(wr);
            }

            //terrainMapImage.copyData(terrainMapWR);
            for (int x = 0; x < dimension1; x++) {
                for (int y = 0; y < dimension2; y++) {
                    for (ImageFile file : files) {
                        double value = file.convert(ccd[y][x].amplitude(), min, max);
                        file.getWritableRaster().setPixel(x, y, getValue(value));
                    }

                }
            }
            for (ImageFile file : files) {
                file.getBufferedImage().setData(file.getWritableRaster());
                ImageIO.write(file.getBufferedImage(), formatName, new File(outputFolder, file.getFileName() + counter + ".png"));
            }
            counter++;
        } catch (IOException e) {
            e.printStackTrace();
        }
    }

    private int[] getValue(double value) {
        return new int[]{(int) (255 * value), (int) (255 * value), (int) (255 * value), (int) (255 * value)};
    }

}

CurveFitting.java

package com.clarke.agnes.experiments;

import java.util.Random;

public class CurveFitting {

    public static final double SAMPLE_SEPARATION = 0.04;
    public static final Random random = new Random();
    public static final double FINENESS = 20.0;
    private static final double WIDE_STEP = Math.PI / FINENESS;
    private static final double MEDIUM_STEP = WIDE_STEP / FINENESS;
    private static final double FINE_STEP = MEDIUM_STEP / FINENESS;
    private static final double FINEST_STEP = FINE_STEP / FINENESS;
    private static final int ITERATIONS = 200 * 200;

    private static  int stepsToSolve = 0;
    private static double errorFit, errorRough, errorRoughHr;
    private static double errorFitMax, errorRoughMax, errorRoughHrMax;
    private static double duration = 0;

    private final double sampleSeparation;

    public CurveFitting(double sampleSeparation) {
        this.sampleSeparation = sampleSeparation;
    }

    public static void main(String[] args) {
        long start = System.nanoTime();
        for (int i = 0; i < ITERATIONS; i++) {
            execute();
        }
        long end = System.nanoTime();
        System.out.printf("Avg Error by curve fitting        : %.5f%s\n", errorFit / ITERATIONS, '%');
        System.out.printf("Avg Error by rough curve sampling : %.5f%s\n", errorRough / ITERATIONS, '%');
        System.out.printf("Avg Error by fine curve sampling  : %.5f%s\n", errorRoughHr / ITERATIONS, '%');
        System.out.println();
        System.out.printf("Max Error by curve fitting        : %.5f%s\n", errorFitMax, '%');
        System.out.printf("Max Error by rough curve sampling : %.5f%s\n", errorRoughMax, '%');
        System.out.printf("Max Error by fine curve sampling  : %.5f%s\n", errorRoughHrMax, '%');
        System.out.println();
        System.out.println("Steps for fine solution           : " + stepsToSolve / ITERATIONS);
        System.out.printf("Millis taken for ccd              : %.1f\n", duration);
        System.out.printf("Millis taken for program          : %.1f\n", (end - start) / 1000000.0);
    }

    private static void execute() {
        double amplitude = random.nextDouble() * 3.0;
        double phase = random.nextDouble() * Math.PI * 2.0;
        Sinusoid unknownWave = new Sinusoid(1.0, amplitude, phase);
        //rough estimate
        double roughEstimateAmp = 0.0;
        for (double xx = 0.0; xx < Math.PI * 2.0; xx += Math.PI / 6.0) {
            double smp = Math.abs(unknownWave.getSample(xx).y);
            roughEstimateAmp = roughEstimateAmp > smp ? roughEstimateAmp : smp;
        }
        double roughEstimateAmpHiRes = 0.0;
        for (double xx = 0.0; xx < Math.PI; xx += Math.PI / 12.0) {
            double smp = Math.abs(unknownWave.getSample(xx).y);
            roughEstimateAmpHiRes = roughEstimateAmpHiRes > smp ? roughEstimateAmpHiRes : smp;
        }


        //get it exactly
        Point[] samples = get2samples(unknownWave);

        Point sample1 = samples[0];
        Point sample2 = samples[1];

        long start = System.nanoTime();
        CurveFitting fitter = new CurveFitting(SAMPLE_SEPARATION);
        Fit bestFit4 = fitter.fitCurve(sample1, sample2);
        long end = System.nanoTime();
        duration += (end - start) / 1000000.0;

        double errorFitLocal = (Math.abs(1.0 - bestFit4.amplitudeSolution / unknownWave.amplitude)) * 100.0;
        double errorRoughLocal = (Math.abs(1.0 - roughEstimateAmp / unknownWave.amplitude)) * 100.0;
        double errorRoughHrLocal = (Math.abs(1.0 - roughEstimateAmpHiRes / unknownWave.amplitude)) * 100.0;
        if (errorFitLocal > 1.0) {
            System.out.printf("Solution3 for phase = %.7f (target=%.7f, %.5f)\n", bestFit4.phaseSolution, unknownWave.phase, (Math.abs(1.0 - bestFit4.phaseSolution/unknownWave.phase)) * 100.0);
            System.out.printf("Solution3 for ampli = %.7f (target=%.7f, %.5f)\n", bestFit4.amplitudeSolution, unknownWave.amplitude, errorFitLocal);
            System.out.printf("RoughEsti for ampli = %.7f (target=%.7f, %.5f)\n", roughEstimateAmp, unknownWave.amplitude, errorRoughLocal);
            System.out.printf("RoughEsHR for ampli = %.7f (target=%.7f, %.5f)\n", roughEstimateAmpHiRes, unknownWave.amplitude, errorRoughHrLocal);
            System.out.println();
        }
        errorFit += errorFitLocal;
        errorRough += errorRoughLocal;
        errorRoughHr += errorRoughHrLocal;
        errorFitMax = errorFitLocal > errorFitMax ? errorFitLocal : errorFitMax;
        errorRoughMax = errorRoughLocal > errorRoughMax ? errorRoughLocal : errorRoughMax;
        errorRoughHrMax = errorRoughHrLocal > errorRoughHrMax ? errorRoughHrLocal : errorRoughHrMax;
    }

    public Fit fitCurve(Point sample1, Point sample2) {
        Fit bestFit = fitCurve(sample1, sample2, 0, Math.PI * 2.0, WIDE_STEP);
        Fit bestFit2 = fitCurve(sample1, sample2, bestFit.xSolution - WIDE_STEP, bestFit.xSolution + (2 * WIDE_STEP), MEDIUM_STEP);
        Fit bestFit3 = fitCurve(sample1, sample2, bestFit2.xSolution - MEDIUM_STEP, bestFit2.xSolution + (2 * MEDIUM_STEP), FINE_STEP);
        return fitCurve(sample1, sample2, bestFit3.xSolution - FINE_STEP, bestFit3.xSolution + (2 * FINE_STEP), FINEST_STEP);
    }

    private static Point[] get2samples(Sinusoid unknownWave) {
        Point[] points = {
                unknownWave.getSample(0.1),
                unknownWave.getSample(0.1 + SAMPLE_SEPARATION),
                unknownWave.getSample(0.1 + SAMPLE_SEPARATION * 2),
                unknownWave.getSample(0.1 + SAMPLE_SEPARATION * 3)
        };
        return get2samples(points);
    }

    public static Point[] get2samples(Point[] points) {
        Point minimumAbs = points[0];
        for (Point point : points) {
            double y = Math.abs(point.y);
            if (y < Math.abs(minimumAbs.y)) {
                minimumAbs = point;
            }
        }
        for (int i = 0; i < points.length - 1; i++) {
            Point p1 = points[i];
            Point p2 = points[i + 1];
            if (minimumAbs.equals(p1) || minimumAbs.equals(p2)) {
                continue;
            }
            return new Point[]{p1, p2};
        }
        throw new RuntimeException("Unreachable statement");
    }

    private Fit fitCurve(Point sample1, Point sample2, double from, double to, double step) {
        //Get lucky with divbyzero
        double yRatio = sample1.y / sample2.y;
        boolean y1Positive = sample1.y >= 0.0;

        Fit bestFit = new Fit(0, 0, Double.MAX_VALUE, 0);

        //y = (amp)sin(x + phase)
        for (double x = from; x < to; x += step) {
            stepsToSolve++;
            double y1 = Math.sin(x * Math.PI * 2.0);
            double y2 = Math.sin((x + sampleSeparation) * Math.PI * 2.0);
            double yRatio2 = y1 / y2;
            boolean y1Positive2 = y1 >= 0.0;
            if (y2 == 0 || sample2.y == 0 || yRatio2 == 0.0) {
                continue;
            }
            double fit = Math.abs(1.0 - yRatio/yRatio2);
            if (y1Positive == y1Positive2 && fit < bestFit.fit ) {
                bestFit.phaseSolution = x - sample1.x;
                bestFit.amplitudeSolution = (sample1.y / y1);
                bestFit.fit = fit;
                bestFit.xSolution = x;
            }
        }
        return bestFit;
    }

    public static class Point {
        public double x;
        public double y;
        public Point(double x, double y) {
            this.x = x;
            this.y = y;
        }
    }

    public static class Sinusoid {
        public double wavelength;
        public double amplitude;
        public double phase;
        public Sinusoid(double frequency, double amplitude, double phase) {
            this.wavelength = frequency;
            this.amplitude = amplitude;
            this.phase = phase;
        }
        private Point getSample(double x) {
            double y = Math.sin((x + phase) * Math.PI * 2.0) * amplitude;
            return new Point(x, y);
        }

    }

    public static class Fit {
        public double phaseSolution = 0;
        public double amplitudeSolution = 0;
        public double fit = 1.0;
        public double xSolution = 0;
        public Fit(double phaseSolution, double amplitudeSolution, double fit, double xSolution) {
            this.phaseSolution = phaseSolution;
            this.amplitudeSolution = amplitudeSolution;
            this.fit = fit;
            this.xSolution = xSolution;
        }
    }

}

ImageFile.java

package com.clarke.agnes.experiments;

import java.awt.image.BufferedImage;
import java.awt.image.WritableRaster;

/**
 * Converts photosite signal to a pixel intensity.
 * Created by Ags on 4/10/2015.
 */
public class ImageFile {

    private String fileName;
    private Converter converter;
    private BufferedImage bufferedImage;
    private WritableRaster writableRaster;

    public ImageFile(String fileName, Converter converter) {
        this.fileName = fileName;
        this.converter = converter;
    }

    public String getFileName() {
        return fileName;
    }

    public double convert(double amplitude, double min, double max) {
        return converter.convert(amplitude, min, max);
    }

    public void setBufferedImage(BufferedImage bufferedImage) {
        this.bufferedImage = bufferedImage;
    }

    public void setWritableRaster(WritableRaster writableRaster) {
        this.writableRaster = writableRaster;
    }

    public BufferedImage getBufferedImage() {
        return bufferedImage;
    }

    public WritableRaster getWritableRaster() {
        return writableRaster;
    }

    public static interface Converter {
        public double convert(double amplitude, double min, double max);
    }
}

TelescopeSimulator.java

package com.clarke.agnes.experiments;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;
import java.util.*;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;

/**
 * Simulates a telescope.
 * Created by Ags on 3/21/2015.
 */
public class TelescopeSimulator {

    private static final int THREADS = 8;

    private static final double RADIANS_TO_DEGREES = 57.295779513082320876798154814105;
    public static final double NM = 1.0 / 1000 / 1000 / 1000;
    public static final double WAVELENGTH_GREEN = 510 * NM;
    public static Random R = new Random(1234567L);
    public static List<Point> wavefront = new ArrayList<>();
    public static List<Point> optic = new ArrayList<>();
    private static int sampleCounter;

    public static void simulate(double focalLength, double aperture, double centralObstructionDiameterRatio, double vaneWidth,
                                double chipWidthArcSeconds, int pixels, int subs, boolean reflect1, boolean reflect2,
                                ImageFile[] files, File outputFolder) throws InterruptedException, IOException {
        final ExecutorService executorService = Executors.newFixedThreadPool(THREADS);
        long begin = System.nanoTime();
        //4fy = x2
        //4fz = x2 + y2
        // 4 * FOCAL_LENGTH * z = x2 + y2
        CcdToImage img = new CcdToImage();
        Point focus = new Point(0, 0, focalLength);

        double arcSecondsPerPixel = chipWidthArcSeconds / pixels;
        double chipWidthMeters = arcSecondsPerPixel / RADIANS_TO_DEGREES * focalLength * pixels / 3600;
        System.out.println("Chip width in mm = " + (chipWidthMeters * 1000.0));

        //generateccd
        Photosite[][] ccd = new Photosite[pixels][pixels];
        for (int i = 0; i < pixels; i++) {
            for (int k = 0; k < pixels; k++) {
                double x = focus.x - (chipWidthMeters / 2.0) + (chipWidthMeters * i / pixels);
                double y = focus.y - (chipWidthMeters / 2.0) + (chipWidthMeters * k / pixels);
                ccd[i][k] = new Photosite(x, y, focus.z);
            }
        }

        //generate resolution finder
        Photosite[] resolutionFinder = new Photosite[(int) (20 / aperture)];
        double metersPerHundredthArcSecond = chipWidthMeters / chipWidthArcSeconds / 100.0;

        for (int i = 0; i < resolutionFinder.length; i++) {
            double x = focus.x + (i * metersPerHundredthArcSecond);
            double y = focus.y;
            resolutionFinder[i] = new Photosite(x, y, focus.z);
        }

        Collection<Callable<Object>> scanners = new ArrayList<>(THREADS);
        for (int i = 0; i < THREADS; i++) {
            scanners.add(new ParallelScanner(THREADS, i, pixels, ccd, reflect1, reflect2));
        }

        double sampleSpacing = aperture / 100.0;
        double halfSampleSpacing = sampleSpacing / 2.0;

        for (int sub = 0; sub < subs; sub++) {
            wavefront = new ArrayList<>();
            optic = new ArrayList<>();
            System.out.println("Image " + sub);
            //generate wavefront and optic
            for (double x = 0 - (aperture / 2.0); x < (aperture / 2.0); x += sampleSpacing) {
                for (double y = 0 - (aperture / 2.0); y < (aperture / 2.0); y += sampleSpacing) {
                    double offsetX = (sampleSpacing * R.nextDouble()) - halfSampleSpacing;
                    double offsetY = (sampleSpacing * R.nextDouble()) - halfSampleSpacing;
                    double xActual = x + offsetX;
                    double yActual = y + offsetY;
                    if (xActual * xActual + yActual * yActual > (aperture / 2.0) * (aperture / 2.0)) {
                        continue;
                    }
                    if (xActual * xActual + yActual * yActual < (aperture * centralObstructionDiameterRatio / 2.0) * (aperture * centralObstructionDiameterRatio / 2.0)) {
                        continue;
                    }
                    if (Math.abs(yActual) < vaneWidth / 2.0) {
                        continue;
                    }
                    if (Math.abs(xActual) < vaneWidth / 2.0) {
                        continue;
                    }
                    Point start = new Point(xActual, yActual, focalLength * 2);
                    double z = (xActual * xActual + yActual * yActual) / (4.0 * focalLength);
                    Point incidence = new Point(xActual, yActual, z);
                    wavefront.add(start);
                    optic.add(incidence);
                }
            }
            scanResolutionFinder(resolutionFinder);
            executorService.invokeAll(scanners);
            reflection(ccd, reflect1, reflect2);
            StringBuilder sb = new StringBuilder();
            sb.append(setupData(focalLength, aperture, centralObstructionDiameterRatio, vaneWidth, chipWidthArcSeconds, pixels, subs));
            int mid = (int) (pixels / 2.0);
            double max = 0, min = Double.MAX_VALUE;
            for (int i = 0; i < pixels; i++) {
                max = max > ccd[i][mid].amplitude() ? max : ccd[i][mid].amplitude();
                min = min < ccd[i][mid].amplitude() ? min : ccd[i][mid].amplitude();
            }
            writeCcdToFile(ccd, sb);
            double minimumRadius = findResolution(resolutionFinder, sb);
            System.out.println("Resolution in arc seconds: " + minimumRadius);
            analysePattern(ccd, mid, minimumRadius, arcSecondsPerPixel, outputFolder);

            FileOutputStream fos = new FileOutputStream(new File(outputFolder, "diffraction.data"));
            fos.write(sb.toString().getBytes());
            fos.close();
            img.noiseToImage(files, ccd, max, min, "PNG", outputFolder);

        }

        long end = System.nanoTime();
        System.out.printf("Took %.3fms.\n", ((end - begin) / 1000000.0));
        executorService.shutdown();
    }

    static final StringBuilder analysis = new StringBuilder(
            "Resolution,Airy Disc Energy,X2 Radii,X3 Radii,X4 Radii,X5Radii,Contrast X2,Contrast X3,Contrast X4,Contrast X5,Contrast 8 arcsec,Contrast 30 arcsec\r\n");
    private static void analysePattern(Photosite[][] ccd, int mid, double minimumRadius, double arcSecondsPerPixel, File outputFolder) throws IOException {
        double airyRadiusInPixels = minimumRadius / arcSecondsPerPixel;
        System.out.println("Airy radius in pixels: " + airyRadiusInPixels);
        if (airyRadiusInPixels < 3) {
            throw new RuntimeException("Pattern is undersampled, airy disc is only " + airyRadiusInPixels + " pixies in radius.");
        }
        double airyDisc = integrate(ccd, mid, airyRadiusInPixels);
        double airyDiscX2 = integrate(ccd, mid, airyRadiusInPixels * 2) - airyDisc;
        double airyDiscX3 = integrate(ccd, mid, airyRadiusInPixels * 3) - airyDisc - airyDiscX2;
        double airyDiscX4 = integrate(ccd, mid, airyRadiusInPixels * 4) - airyDisc - airyDiscX2 - airyDiscX3;
        double airyDiscX5 = integrate(ccd, mid, airyRadiusInPixels * 5) - airyDisc - airyDiscX2 - airyDiscX3 - airyDiscX4;
        double eightASDisc = integrate(ccd, mid, (int) (8.0 / arcSecondsPerPixel));
        double thirtyAsDisc = integrate(ccd, mid, (int) (30.0 / arcSecondsPerPixel));
        System.out.println("Energy at Airy disc radii:");

        double ad = airyDisc / airyDisc;
        double adx2 = airyDiscX2 / airyDisc;
        double adx3 = airyDiscX3 / airyDisc;
        double adx4 = airyDiscX4 / airyDisc;
        double adx5 = airyDiscX5 / airyDisc;
        System.out.println("\tDisc Radius x1: " + ad);
        System.out.println("\tDisc Radius x2: " + adx2);
        System.out.println("\tDisc Radius x3: " + adx3);
        System.out.println("\tDisc Radius x4: " + adx4);
        System.out.println("\tDisc Radius x5: " + adx5);
        System.out.println("Airy disc to surroundings energy ratio (contrast):");

        double contrastx2 = airyDisc / (airyDisc + airyDiscX2);
        double contrastx3 = airyDisc / (airyDisc + airyDiscX2 + airyDiscX3);
        double contrastx4 = airyDisc / (airyDisc + airyDiscX2 + airyDiscX3 + airyDiscX4);
        double contrastx5 = airyDisc / (airyDisc + airyDiscX2 + airyDiscX3 + airyDiscX4 + airyDiscX5);
        double eightAsContrast = airyDisc / eightASDisc;
        double thirtyAsContrast = airyDisc / thirtyAsDisc;
        System.out.println("\tAiry disc v Airy disc x2: " + contrastx2);
        System.out.println("\tAiry disc v Airy disc x3: " + contrastx3);
        System.out.println("\tAiry disc v Airy disc x4: " + contrastx4);
        System.out.println("\tAiry disc v Airy disc x5: " + contrastx5);
        System.out.println("\tAiry disc v 4 arc second radius disc: " + eightAsContrast);
        System.out.println("\tAiry disc v Jupiter disc: " + thirtyAsContrast);

        analysis.append(minimumRadius).append(',');
        analysis.append(ad).append(',').append(adx2).append(',').append(adx3).append(',').append(adx4).append(',').append(adx5).append(',');
        analysis.append(contrastx2).append(',').append(contrastx3).append(',').append(contrastx4).append(',').append(contrastx5).append(',');
        analysis.append(eightAsContrast).append(',').append(thirtyAsContrast).append('\r').append('\n');

        FileOutputStream fos = new FileOutputStream(new File(outputFolder, "analysis.csv"));
        fos.write(analysis.toString().getBytes());
        fos.close();
    }

    private static double integrate(Photosite[][] ccd, int mid, double radius) {
        double radiusSquared = radius * radius;
        double result = 0.0;
        for (int i = mid - (int)radius; i < mid + radius; i++) {
            for (int k = mid - (int)radius; k < mid + radius; k++) {
                double distanceSquared = (i - mid) * (i - mid) * (k - mid) * (k - mid);
                if (distanceSquared > radiusSquared) {
                    continue;
                }
                result += ccd[i][k].amplitude() * ccd[i][k].amplitude();
            }
        }
        return result;
    }

    private static void writeCcdToFile(Photosite[][] ccd, StringBuilder sb) {
        for (Photosite[] photosites : ccd) {
            for (Photosite photosite : photosites) {
                boolean first = true;
                for (double sample : photosite.samples) {
                    if (first) {
                        first = false;
                    }else {
                        sb.append('|');
                    }
                    sb.append(sample);
                }
                sb.append('\n');
            }
        }
    }

    private static double findResolution(Photosite[] resolutionFinder, StringBuilder sb) {
        double distanceFromCenter = 0.0;
        double minimumRadius = 0.0;
        boolean foundUptick = false;
        double previousAmplitude = resolutionFinder[0].amplitude();
        for (Photosite photosite : resolutionFinder) {
            sb.append(String.format("%.2f\t%s\n", distanceFromCenter, photosite.amplitude()));
            double currentAmplitude = photosite.amplitude();
            if (currentAmplitude > previousAmplitude) {
                foundUptick = true;
            }
            if (!foundUptick) {
                previousAmplitude = currentAmplitude;
                minimumRadius = distanceFromCenter;
            }
            distanceFromCenter += 0.01;
        }
        if (!foundUptick) {
            throw new RuntimeException("No uptick found in resolutionFinder, check the radius of resolutionFinder.");
        }
        return minimumRadius;
    }

    private static void scanResolutionFinder(Photosite[] resolutionFinder) {
        for (Photosite p : resolutionFinder) {
            for (int j = 0; j < wavefront.size(); j++) {
                Point start = wavefront.get(j);
                Point incidence = optic.get(j);
                double fromPlane = start.distance(incidence);
                //Illuminate ccd
                double distance = incidence.distance(p) + fromPlane;
                double phase = distance % WAVELENGTH_GREEN;
                int counter = 0;
                for (double t = 0; t < 0.15; t += CurveFitting.SAMPLE_SEPARATION) {
                    double phaseSample = (((WAVELENGTH_GREEN * t) + phase) % WAVELENGTH_GREEN) / WAVELENGTH_GREEN;
                    double amplitude = Math.sin(phaseSample * 2 * Math.PI);
                    p.samples[counter++] += amplitude;
                }
            }
        }
    }

    private static void reflection(Photosite[][] ccd, boolean reflect1, boolean reflect2) {
        int mid = ccd.length / 2;
        if (reflect1) {
            for (int i = mid; i < ccd.length; i++) {
                int reflected = ccd.length - i;
                ccd[reflected] = Arrays.copyOf(ccd[i], ccd.length);
            }
        }
        if (reflect2) {
            for (int i = 0; i < ccd.length; i++) {
                for (int k = mid; k < ccd.length; k++) {
                    int reflected = ccd.length - k;
                    ccd[i][reflected] = ccd[i][k];
                }
            }
        }
    }

    private static String setupData(double focalLength, double aperture, double centralObstructionDiameterRatio,
                                    double vaneWidth, double chipWidthArcSeconds, int pixels, int subs) {
        return "focalLength:" + focalLength +
                "\naperture:" + aperture +
                "\ncentralObstructionDiameterRatio:" + centralObstructionDiameterRatio +
                "\nvaneWidth:" + vaneWidth +
                "\nchipWidthArcSeconds:" + chipWidthArcSeconds +
                "\npixels:" + pixels +
                "\nsubs:" + subs + "\n-----------------------------------------\n";
    }

    public static class Point {
        public final double x;
        public final double y;
        public final double z;
        private Point distanceFromPoint = null;
        private double distanceFromDistance = 0.0;

        public Point(double x, double y, double z) {
            this.x = x;
            this.y = y;
            this.z = z;
        }

        public double distance(Point p) {
            if (p.equals(distanceFromPoint)) {
                return distanceFromDistance;
            }
            distanceFromPoint = p;
            distanceFromDistance = Math.sqrt(((p.x - x) * (p.x - x)) + ((p.y - y) * (p.y - y)) + ((p.z - z) * (p.z - z)));
            return distanceFromDistance;
        }
    }

    public static class Photosite extends Point {

        public double[] samples = new double[4];
        private static final CurveFitting fitting = new CurveFitting(CurveFitting.SAMPLE_SEPARATION);

        public Photosite(double x, double y, double z) {
            super(x, y, z);
        }

        public double amplitude() {
            CurveFitting.Point[] points = new CurveFitting.Point[]{
                    new CurveFitting.Point(0.0, samples[0]),
                    new CurveFitting.Point(CurveFitting.SAMPLE_SEPARATION, samples[1]),
                    new CurveFitting.Point(CurveFitting.SAMPLE_SEPARATION * 2.0, samples[2]),
                    new CurveFitting.Point(CurveFitting.SAMPLE_SEPARATION * 3.0, samples[3])
            };
            points = CurveFitting.get2samples(points);
            CurveFitting.Fit fit = fitting.fitCurve(points[0], points[1]);
            return fit.amplitudeSolution;
        }
    }

    public static class ParallelScanner implements Callable<Object> {

        private final int step;
        private final int offset;
        private final int pixels;
        private final Photosite[][] ccd;
        private final boolean reflect1;
        private final boolean reflect2;

        public ParallelScanner(int step, int offset, int pixels, Photosite[][] ccd, boolean reflect1, boolean reflect2) {
            this.step = step;
            this.offset = offset;
            this.pixels = pixels;
            this.ccd = ccd;
            this.reflect1 = reflect1;
            this.reflect2 = reflect2;
        }

        @Override
        public Object call() throws Exception {
            for (int i = offset; i < pixels; i += step) {
                int mid = (int) (pixels / 2.0);
                if (offset == 0) System.out.print("Row " + sampleCounter++ + "\r");
                if (i < mid && reflect1) {
                    continue;
                }
                for (int k = 0; k < pixels; k++) {
                    if (k < mid && reflect2) {
                        continue;
                    }
                    Photosite p = ccd[i][k];
                    for (int j = 0; j < wavefront.size(); j++) {
                        Point start = wavefront.get(j);
                        Point incidence = optic.get(j);
                        double fromPlane = start.distance(incidence);
                        //Illuminate ccd
                        double distance = incidence.distance(p) + fromPlane;
                        double phase = distance % WAVELENGTH_GREEN;
                        int counter = 0;
                        for (double t = 0; t < 0.15; t += CurveFitting.SAMPLE_SEPARATION) {
                            double phaseSample = (((WAVELENGTH_GREEN * t) + phase) % WAVELENGTH_GREEN) / WAVELENGTH_GREEN;
                            double amplitude = Math.sin(phaseSample * 2 * Math.PI);
                            p.samples[counter++] += amplitude;
                        }
                    }
                }
            }
            return null;
        }
    }

}

TelescopeSimulatorMain.java

package com.clarke.agnes.experiments;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.IOException;

/**
 * Runs the telescope simulator.
 * Created by Ags on 4/10/2015.
 */
public class TelescopeSimulatorMain {

    public static void main(String[] args) throws InterruptedException, IOException {
        File outputFolder = new File("C:\\Users\\JoeBloggs\\Desktop\\diffraction");

        ImageFile[] files = new ImageFile[]{
                new ImageFile("amplitude", new ImageFile.Converter() {
                    @Override
                    public double convert(double amplitude, double min, double max) {
                        return amplitude / max;
                    }
                }),
                new ImageFile("sqrt", new ImageFile.Converter() {
                    @Override
                    public double convert(double amplitude, double min, double max) {
                        return Math.sqrt(Math.sqrt(amplitude / max));
                    }
                }),
                new ImageFile("eye", new ImageFile.Converter() {
                    double[] level1 = {0.2, 0.5}; //take 0.1 brightness to 0.5
                    double[] level2 = {0.75, 0.6};
                    double[] level3 = {1.0, 1.0};
                    @Override
                    public double convert(double amplitude, double min, double max) {
                        double intensity = amplitude / max;
                        if (intensity <= level1[0]) {
                            intensity = intensity / level1[0] * level1[1];
                        } else if (intensity <= level2[0]) {
                            intensity = level1[1] + ((intensity - level1[0]) / (level2[0] - level1[0]) * (level2[1] - level1[1]));
                        } else {
                            intensity = level2[1] + ((intensity - level2[0]) / (level3[0] - level2[0]) * (level3[1] - level2[1]));
                        }
                        return intensity;
                    }
                })
        };
//        File outputFolderMn190 = new File(outputFolder, "MN190"); outputFolderMn190.mkdir();
//        TelescopeSimulator.simulate(1.0, 0.190, 0.2684, 0.0, 60.0, 500, 16, true, true, files, outputFolderMn190);
//
//        File outputFolderSE4 = new File(outputFolder, "SE4"); outputFolderSE4.mkdir();
//        TelescopeSimulator.simulate(1.324, 0.102, 0.3302, 0.0, 60.0, 500, 16, true, true, files, outputFolderSE4);
//
//        File outputFolderAPM180 = new File(outputFolder, "APM180"); outputFolderAPM180.mkdir();
//        TelescopeSimulator.simulate(1.260, 0.180, 0.0 , 0.0, 60.0, 500, 16, true, true, files, outputFolderAPM180);

//        File outputFolderC11 = new File(outputFolder, "C11"); outputFolderC11.mkdir();
//        TelescopeSimulator.simulate(2.8, 0.280, 0.34, 0.0, 60.0, 500, 16 , true, true, files, outputFolderC11);
        
//        File outputFolderXXXXXX = new File(outputFolder, "XXXXXX"); outputFolderXXXXXX.mkdir();
//        TelescopeSimulator.simulate(1.260, 0.180, 0.0 , 0.0, 60.0, 500, 16, true, true, files, outputFolderXXXXXX);

        File outputFolderSW200PF6 = new File(outputFolder, "SW200PF6"); outputFolderSW200PF6.mkdir();
        TelescopeSimulator.simulate(1.200, 0.200, 0.25 , 0.0005, 60.0, 500, 12, true, true, files, outputFolderSW200PF6);

        File outputFolderCFGSO200F4 = new File(outputFolder, "CFGSO200F4"); outputFolderCFGSO200F4.mkdir();
        TelescopeSimulator.simulate(1.200, 0.200, 0.35 , 0.0005, 60.0, 500, 12, true, true, files, outputFolderCFGSO200F4);

        File outputFolderGSO200F5 = new File(outputFolder, "GSO200F5"); outputFolderGSO200F5.mkdir();
        TelescopeSimulator.simulate(1.260, 0.180, 0.315 , 0.001, 60.0, 500, 12, true, true, files, outputFolderGSO200F5);

        File outputFolderST80 = new File(outputFolder, "ST80"); outputFolderST80.mkdir();
        TelescopeSimulator.simulate(0.400, 0.080, 0.0 , 0.0, 60.0, 500, 12, true, true, files, outputFolderST80);

        File outputFolderAPMTMB105_650 = new File(outputFolder, "APMTMB105_650"); outputFolderAPMTMB105_650.mkdir();
        TelescopeSimulator.simulate(0.650, 0.105, 0.0 , 0.0, 60.0, 500, 12, true, true, files, outputFolderAPMTMB105_650);
    }

}