Simulated Dawes Limit images

The images below simulate diffraction images from telescopes of 0, 10, 20, and 35% central obstructions, and the appearance of a binary star at the Dawes limit separation with a magnitude difference of 0, 1, 2, and 3mv.

They are generated with an algorithm first proposed by Dennis Di Cicco around 20+ years ago in a Sky and Telescope article. It's a rather brute force approach that divides an objective into a series of point sources on a sphere centered on the image plane and optical axis of the objective. In the case of this simulation, I chose a radius of 200 points, giving 125,663 points for the objective.

Initially, the dimensions of the aperture and focal length are in centimeters. This is converted to units of 500 nanometers, which is the size of a wavelength of light, close to the maximum visual acuity of the human eye. The heart of this algorithm is that for each of the 320,000 points on the image plane, we calculate the distance to each point on the objective, find out how much it differs from a whole wavelength, and take the sine of this difference and add it to the amplitude of that point. Basic physical optics. But that works out to a bit over 40 billion calculations for each aperture. (It took 4.5 hours to run on my machine.) The image plane is about 13 Dawes' limits in linear size. This works out to approximately an image plane of 74 microns square for my 20cm f/10. Tiny. To see the single star patterns below, you'd have to have a power of over 2000.

Once this is done, we calculate the actual intensity of the light at that point on the focal plane. This is simply the amplitude squared. This quantity is multiplied by a scale factor that maps the brightest pixel to 255. This is then converted to a .jpg image by Java's BufferedImage and ImageIO methods. This final image represents what an ideal single star looks like at very high magnification.

Next, the binary images are formed. To do this the image is then combined with a copy of itself spaced a Dawes' limit from the other. The copies are 0, 1, 2, and 3 magnitudes fainter than the "primary". Note how easily stars of equal magnitudes can be split. And how it's well nigh impossible to split those where the companion is 3 magnitudes fainter than its primary. Going over my own observations, this is correct. I've split some very faint pairs that were only 0.5" apart with my C-8, but had little success with doubles that are over twice a Dawes' limit apart but with a 2-3mv difference.

Note that this simulation is as good as it gets, optically. There are No aberations, the optics are perfectly collimated, the atmosphere is Pickering 10, and the star is monochromatic. It'll never get this good in real life.

I've included a trace of the intensity superimposed on the single star images. It is scaled such that the most intense pixels are plotted a the top of the images, and the least intense pixels are at the bottom.

The Java code that created the images follows the pictures. Feel free to improve it!

Simulated diffraction images from telescopes of 0, 10, 20, and 35% central obstructions.

Unobstructed image.

Magnitude difference of 0.	Magnitude difference of 1.
Magnitude difference of 2.	Magnitude difference of 3.

10% obstruction.

Magnitude difference of 0.	Magnitude difference of 1.
Magnitude difference of 2.	Magnitude difference of 3.

20% obstruction.

Magnitude difference of 0.	Magnitude difference of 1.
Magnitude difference of 2.	Magnitude difference of 3.

35% obstruction.

Magnitude difference of 0.	Magnitude difference of 1.
Magnitude difference of 2.	Magnitude difference of 3.

Here's the Java code that created the above images:

import java.awt.Color;
import java.awt.image.BufferedImage;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
// import java.io.FileWriter;
import java.io.IOException;

import javax.imageio.ImageIO;

public class Diffraction {
 
  // Calculate the diffraction pattern for a series of telescopes,
  // with no abberations and 0 to 35% central obstructions.
  
  public void calc(double aperture,
                   double fRatio) {
    aperture *= 2e4;  // The aperture, in wavelength units.
    
    BufferedImage bid[] = new BufferedImage[24];
    int bidCt = 0;
    
    // The binary image is based on combining two 400 pixel images.
    int nrPxls = 400,                           
        halfPix = nrPxls / 2,
        dPix = (int) (nrPxls * (Math.sqrt(2))), // Dimension of the diffraction pattern image.
        hdPix = dPix / 2,
        hp2 = halfPix * halfPix;
   
    // The image field is 6 wavelengths across. (75 microns for a 20cm f/10)
    double apIncr = (aperture / (double) nrPxls),
           f = aperture * fRatio,
           f2 = f * f,
           fpIncr = aperture * fRatio * 6 / (206264.81 * nrPxls);
    
    double image[][] = new double[dPix + 1][dPix + 1 ];
 
    for (int s = 0; s < 4; s++) {
      double obs = s;
      if (s == 3) { obs = 3.5; } // The obstruction of an SCT.

      double rs = obs * 0.1 * aperture / 2; // Radius of the secondary, in
                                            // wavelengths.
      // Traverse the aperture.
      for (int x = -nrPxls; x < nrPxls; x++) {
        double xa = x * apIncr;
        int maxY = (int) (Math.sqrt(hp2 - (x * x)));
        for (int y = -maxY; y < maxY; y++) {
          double ya = (double) y * apIncr;
          // If the point is blocked by the secondary, ignore it.
          double rp = Math.sqrt((xa * xa) + (ya * ya));
          if (rp < rs) { continue; } 

          double z = Math.sqrt(f2 - ((xa * xa) + (ya * ya)));

          // And the image plane.
          for (int i = -(hdPix - 1); i < hdPix; i++) {
            double xi = (double) i * fpIncr;
            for (int j = -(hdPix - 1); j < hdPix; j++) {
              double yi = (double) j * fpIncr;
              double dx = xa - xi,
                     dy = ya - yi,
                     r = Math.sqrt((dx * dx) + (dy * dy) + (z * z));
              int fullWaves = (int) r;
              double fractionalWaves = ((r - (double) fullWaves) * 2) - 1;
              double amplitude = Math.sin((1 - fractionalWaves) * Math.PI / 2);
              image[hdPix + i][hdPix + j] += amplitude;
            }
          }
        }
      }

      // Create a jpeg image of the diffraction pattern just calculated.
      
      // What is the brightest point in the image?
      double maxAmp = 0;
      for (int x = 0; x < dPix; x++) {
        for (int y = 0; y < dPix; y++) {
          if (image[x][y] > maxAmp) {
            maxAmp = image[x][y];            
          }
        }
      }

      // Create 256 grey tones, black to white.
      int[] rgb = new int[256];
      for (int cIndex = 0; cIndex < 256; cIndex++) {
        Color c = new Color(cIndex, cIndex, cIndex);
        rgb[cIndex] = c.getRGB();
      }
      
      bid[bidCt] = new BufferedImage(dPix, dPix, BufferedImage.TYPE_INT_RGB);

      int value;                   // The intensity of a given pixel
      double k = 255 / (maxAmp * maxAmp);
      for (int x = 0; x < dPix; x++) {
        for (int y = 0; y < dPix; y++) {
          
          // Intensity is amplitude squared.
          value = (int) ((image[x][y] * image[x][y]) * k);
          
          // Make a white star on a black background, not the reverse.
          value = 255 - value;
          bid[bidCt].setRGB(x, y, rgb[value]);
        }
        
        // Make a white pixel for the graph of the intensity. An intensity of 0 should 
        // plot to the bottom of the image, one of 255 to the top.
        value = (int) (((image[x][hdPix] * image[x][hdPix]) * k) * ((double) dPix / 255));
        bid[bidCt].setRGB(x, value, rgb[255]);
      }

      try {
        ImageIO.write(bid[bidCt], "JPEG", new File("/work/glxy/tmp/SingleStar_" + s + ".jpg"));
      } catch (IOException ie) {
        ie.printStackTrace();
      }

      bidCt++;
      
      // Now combine the image with it's fainter version a Dawes' limit away to simulate
      // the appearance of a double star at the limit of the telescopes resolution. We
      // make 4 images, for magnitude differences of 0 - 3 magnitudes difference.
      
      double dawes = 11.43 * 2e4 / aperture,     // Dawes limit in wavelengths.
      mv = Math.pow(10, 0.4);                    // One magnitude of intensity difference.
      int delta = (dPix - nrPxls) / 2;           // Single star image is larger than double image.
      int offset = (int) (dawes * halfPix / 3);  // Dawes Offset for each diffraction pattern.
      int xMax = nrPxls + offset;
      int yOffset = hdPix - halfPix;
           
      // Calculate images for stars of equal magnitude, and drop a magnitude with each iteration.
      for (int i = 0; i < 4; i++) {
        double dImg[][] = new double[xMax + 1][nrPxls];
        double faintFactor = Math.pow(mv, (double) i);
        maxAmp = 0;        
        for (int x = 0; x < xMax; x++) {
          for (int y = 0; y < nrPxls; y++) {          
            dImg[x][y] += (image[x + delta][y + yOffset]) / faintFactor +
                           image[x + delta + offset][y + yOffset];
  
            if (dImg[x][y] > maxAmp) { maxAmp = dImg[x][y]; }
          }
        }
        
        bid[bidCt] = new BufferedImage((xMax + 1), nrPxls, BufferedImage.TYPE_INT_RGB);
        k = 255 / (maxAmp * maxAmp);
        for (int x = 0; x < xMax; x++) {
          for (int y = 0; y < nrPxls; y++) {
 
            // Intensity is amplitude squared.
            value = (int) ((dImg[x][y] * dImg[x][y]) * k);
            
            // Make a white star on a black background, not the reverse.
            value = 255 - value;
            bid[bidCt].setRGB(x, y, rgb[value]);
          }
        }

        try {
          ImageIO.write(bid[bidCt], "JPEG", new File("/work/glxy/tmp/double_" + i + "_" + s + ".jpg"));
        } catch (IOException ie) {
          ie.printStackTrace();
        }
        
        bidCt++;
      }
    }
  }
  
  static public void main(String args[]) {
    
    long startTime = System.currentTimeMillis();
    
    double aperture = 20,
           fRatio = 10;
    
    if (args.length == 2) {
       aperture = Double.parseDouble(args[0]);
       fRatio = Double.parseDouble(args[1]);
    }
    
    Diffraction d = new Diffraction();
    d.calc(aperture, fRatio);
    
    long endTime   = System.currentTimeMillis();
    int totalTime = (int) (endTime - startTime);
    int min = totalTime / 60000;
    int sec = (totalTime - (min * 60000)) / 1000;
    
    System.out.println("The run took " + min + " minutes, " + sec + " seconds.");
  }
}