Straight Edge Diffraction Pattern from a Point Source

A nice way to find the angular diameters of stars is by Lunar occultations of stars. Since ordinary telescopes only show stars as point objects, this cannot be directly measured by observation. HST and such huge observatories with precise optics alone can resolve such details.

What limits the resolution of a telescope and thus our capability to see the disks of stars is diffraction, but the very same diffraction can help us find out the angular diameter of stars! When the moon passes across a star, the limb of the moon can act like a straight edge (it is straight in comparison to the star’s angular dimensions) and produce a diffraction pattern of the star’s light.

Here is a beautiful article outlining the whole procedure, of finding angular diameters of stars using this technique. Dr. Shylaja at the planetarium wanted me to write a small program to convert such an intensity pattern into the angular diameter of the star. I liked the idea and have started working on it. I’m assuming the formula for the intensity distribution, as he does in the article.

I wrote today, a program to numerically evaluate the intensity pattern from a point source, as explained in the article. I found out that the GNU Scientific Library (GSL) does integrations in a very elegant manner. You needn’t rewrite Simpson’s Rule! Here’s a code snippet that evaluates the required intensity distribution function:

gsl_integration_workspace *w = gsl_integration_workspace_alloc (1000);
double result, error;
double params[2];
double cosint, sinint;
double coserr, sinerr;

params[1] = M_PI / ( screenDist * lambda );
gsl_function F;
F.function = &straightedge_integrand;
F.params = (void *)params;

// TODO: Decide lower limit more intelligently
params[0] = 0;
gsl_integration_qags( &F, -200, -100, 0, 1e-1,
1000, w, &cosint, &coserr );
fprintf(stderr, "Cosine integral from -1000 to -100 = %e\n", cosint);

double b = sqrt( lambda * screenDist );

gsl_integration_qags( &F, -20 * b, x, 0, 1e-7, // TODO: Set appropriately
1000, w, &cosint, &coserr );
params[0] = M_PI / 2;
gsl_integration_qags( &F, -20 * b, x, 0, 1e-7, // TODO: Set appropriately
1000, w, &sinint, &sinerr );

// fprintf(stderr, "cosint = %e, sinint = %e\n", cosint, sinint);
double intensityFactor = cosint * cosint + sinint * sinint;
double relerr = ( ( cosint != 0 && sinint != 0 ) ? 2 * ( coserr / cosint + sinerr / sinint ) : 1.0 );
(*f)[ x ] = intensityFactor;
return intensityFactor;

Somehow, the integral from -inf to x failed to evaluate, probably because the transform that GSL does to evaluate semi-infinite integrals of this form retained the oscillatory behaviour of the sine and cosine. I need to look into this sometime. Instead, I neglected the portion from -inf to -20*b, which is almost zero (intensity deep inside the geometric shadow region).

My program also stores all the function values it keeps in a map. I preferred a hash, but STL has no hashing function for doubles. Probably I should convert double into string using atof and then hash against that, but I don’t know whether the overhead of atof will be acceptable. I guess it will. This kind of storage will, IMO, help immensely when the convolution is performed between the intensity distribution of the star (“1-D kernel”) and the point source diffraction pattern.

I still don’t understand why it should be convolution rather than finding the correlation, i.e. Why is the mirroring of the kernel required? Anyway, both are the same in this case.

Here are today’s results, which have come after some messing around with the STL containers all night: