Java code for Bessel beam simulation


package psf;

import static java.lang.Math.PI;
import static java.lang.Math.asin;
import static java.lang.Math.cos;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
import static java.lang.String.format;
import static math.ApproxBessel.PolApprox_J;
import gnu.trove.TDoubleArrayList;

import java.io.DataOutputStream;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;

public class Psf
{

/**
* TDoubleArrayList test = new TDoubleArrayList(); { test.getQuick(offset); }/
*
* @throws IOException
**/

public static void main(String[] args) throws IOException
{
// Size of computed volume, in pixel.
// To calculate the XY plane, set Z=1
// To calculate the YZ plane, set X=1
// The default center (focus of the beam) is at (0, 0, 0)
// To change the default center, change the variables Xc, Yc, and Zc of doPSF
int Xsize = 100;
int Ysize = 100;
int Zsize = 1;

// Parameters defining the focusing optics
// NA should not be larger than RefIndex, defined in doPsf
// NAin defines the NA of the inner radius of the annulus. Set to 0 to simulate a Gaussian beam
// Alternatively, epsilon defines the ratio between NAin and NA (always <1)
double NA = 0.97;
double NAin = 0;
double epsilon =0;

// To compute more than 1 PSF change totN
int i;
int totN=1;

// Scaling parameters
// wavelength is the wavelength (nm is ok, but can be also um. Needs to be consisten with the scaling parameters defined below)
// Scale_x, scale_y, and scale_z are the scale along the X, Y, and Z directions, respectively (same units as wavelength)
double wavelength = 488;
double scale_x = 10;
double scale_y = 10;
double scale_z = 10;

for(i=0; i<totN; i++){
// Loop used to compute multiple PSF, e. g. at various NAs
//System.out.format("Computing file %d of %d\n", i + 1, totN);
psf.Psf.doPsf(Xsize, Ysize, Zsize, wavelength, scale_x, scale_y, scale_z, NA, epsilon);
NA+=0.01;
}

}

public static void doPsf( int Xsize,
int Ysize,
int Zsize,
double wavelength,
double scale_x,
double scale_y,
double scale_z,
double NA,
double epsilon) throws FileNotFoundException,
IOException
{
{
// Steps to divide the integral into
final int IntSize = 256;

// Ref. index of medium
final double RefIndex = 1;

// PSF center (0,0,0)=volume center
double Xc = 0;
double Yc = 0;
double Zc = 0;

// Wave vector
double kw = 2 * PI / wavelength;

// Max. half-angle defined by NA
double MaxAlpha = java.lang.Math.asin(NA / RefIndex);

// Transversal distance
double r2;

// Axial distance
double z2;

// Angle
double phi;

int i, j, k, l, m;

// Variables for computing elapsed time
// clock_t start;
// clock_t end;
// double deltat;

double arg_Int_cos;
double arg_Int_sin;
double arg_Int_cos_cos;
double arg_Int_sin_cos;
double temp_int_common_R;
double temp_int_common_I;
double temp_int_i0_R;
double temp_int_i0_I;
double temp_int_i1_R;
double temp_int_i1_I;
double temp_int_i2_R;
double temp_int_i2_I;
double temp_int_theta;
// double[] result_array[];

// Starting angle, defined as a fraction of the NA
double AlphaStart = java.lang.Math.asin(epsilon * NA / RefIndex);

// Theta increment for integral calculation
double theta_step = (MaxAlpha - AlphaStart) / IntSize;

// Caching
TDoubleArrayList TheR = new TDoubleArrayList();
TDoubleArrayList J0n = new TDoubleArrayList();
TDoubleArrayList J1n = new TDoubleArrayList();
TDoubleArrayList J2n = new TDoubleArrayList();
int numOfR = 0;
boolean found = false;

System.out.format("Using lambda=%d, NA=%2.2f, e=%2.2f\n",
(int) wavelength,
NA,
epsilon);

numOfR = cacheJ(Xsize,
Ysize,
scale_x,
scale_y,
IntSize,
Xc,
Yc,
kw,
AlphaStart,
theta_step,
TheR,
J0n,
J1n,
J2n,
numOfR,
found);

// end = clock();

// deltat = (end-start)/((double)CLOCKS_PER_SEC);

String I_fileName = format( "%3d_3D_I_%2.2f_%2.2f.raw",
(int) wavelength,
(float) NA,
(float) epsilon);
String Ex_fileName = format("%3d_3D_Ex_%2.2f_%2.2f.raw",
(int) wavelength,
(float) NA,
(float) epsilon);
String Ey_fileName = format("%3d_3D_Ey_%2.2f_%2.2f.raw",
(int) wavelength,
(float) NA,
(float) epsilon);
String Ez_fileName = format("%3d_3D_Ez_%2.2f_%2.2f.raw",
(int) wavelength,
(float) NA,
(float) epsilon);

DataOutputStream I_fos = new DataOutputStream(new FileOutputStream( I_fileName,
true));
DataOutputStream Ex_fos = new DataOutputStream(new FileOutputStream(Ex_fileName,
true));
DataOutputStream Ey_fos = new DataOutputStream(new FileOutputStream(Ey_fileName,
true));
DataOutputStream Ez_fos = new DataOutputStream(new FileOutputStream(Ez_fileName,
true));
// Loops to compute the volume
for (k = 0; k < Zsize; k++)
{

for (j = 0; j < Ysize; j++)
{

for (i = 0; i < Xsize; i++)
{

// Computing the position
r2 = sqrt(((i - Xsize / 2 + Xc)* scale_x) * ((i - Xsize / 2 + Xc)* scale_x)
+ ((j - Ysize / 2 + Yc)* scale_y)
* ((j - Ysize / 2 + Yc)* scale_y));

z2 = (k - Zsize / 2 + Zc) * scale_z;

if (r2 == 0)
{
// Point is on the axis, set phi to 0
phi = 0;
}
else
{
// Point is not on the axis, calculate phi
phi = asin((j - Ysize / 2 + Yc) * scale_y / r2);
}

temp_int_i0_R = 0;
temp_int_i0_I = 0;
temp_int_i1_R = 0;
temp_int_i1_I = 0;
temp_int_i2_R = 0;
temp_int_i2_I = 0;

// Choosing the index m according to r2
// The loop scans TheR, and stops when encounters the right r2.
for (m = 0; m < numOfR; m++)
{
if ((r2 - TheR.getQuick(m)) == 0)
{
// Debug
// System.out.format("Debug: difference is %f\n", (r2-TheR[m]));
break;
}

}

// Debug
// System.out.format("Debug: r2 was %f, chosen %f\n", r2, TheR[m]);

{
for (l = 0; l < IntSize; l++)
{
arg_Int_cos = cos(l * theta_step + AlphaStart);
arg_Int_sin = sin(l * theta_step + AlphaStart);

arg_Int_cos_cos = cos(kw * z2 * arg_Int_cos);
arg_Int_sin_cos = sin(kw * z2 * arg_Int_cos);

// Computing integrals
// See Foreman, M. R., & Török, P. (2011). Computational methods in vectorial imaging. Journal of Modern Optics, 58(5-6), 339–364. doi:10.1080/09500340.2010.525668
// for details

// I(theta). Set to 1 for uniform illumination.
temp_int_theta = 1;//exp(-arg_Int_sin*arg_Int_sin/0.04);
temp_int_common_R = temp_int_theta * theta_step * sqrt(arg_Int_cos)
* arg_Int_cos_cos
* arg_Int_sin;

temp_int_common_I = temp_int_theta * theta_step * sqrt(arg_Int_cos)
* arg_Int_sin_cos
* arg_Int_sin;

temp_int_i0_R += temp_int_common_R * J0n.getQuick(m * IntSize
+ l)
* (1 + arg_Int_cos);

temp_int_i0_I += temp_int_common_I * J0n.getQuick(m * IntSize
+ l)
* (1 + arg_Int_cos);

temp_int_i1_R += temp_int_common_R * J1n.getQuick(m * IntSize
+ l)
* arg_Int_sin;

temp_int_i1_I += temp_int_common_I * J1n.getQuick(m * IntSize
+ l)
* arg_Int_sin;

temp_int_i2_R += temp_int_common_R * J2n.getQuick(m * IntSize
+ l)
* (arg_Int_cos - 1);

temp_int_i2_I += temp_int_common_I * J2n.getQuick(m * IntSize
+ l)
* (arg_Int_cos - 1);
}
}

writeFields(wavelength,
r2,
phi,
temp_int_i0_R,
temp_int_i0_I,
temp_int_i1_R,
temp_int_i1_I,
temp_int_i2_R,
temp_int_i2_I,
I_fos,
Ex_fos,
Ey_fos,
Ez_fos);
}

}

System.out.format("Computing plane %d of %d\n", k + 1, Zsize);

}

I_fos.close();
Ex_fos.close();
Ey_fos.close();
Ez_fos.close();
}
}

private static void writeFields(double wavelength,
double r2,
double phi,
double temp_int_i0_R,
double temp_int_i0_I,
double temp_int_i1_R,
double temp_int_i1_I,
double temp_int_i2_R,
double temp_int_i2_I,
DataOutputStream I_fos,
DataOutputStream Ex_fos,
DataOutputStream Ey_fos,
DataOutputStream Ez_fos) throws IOException
{
double cosphi;
double cos2phi;
double sin2phi;
double Ex_R;
double Ex_I;
double Ey_R;
double Ey_I;
double Ez_R;
double Ez_I;
double I_;
double Ex;
double Ey;
double Ez;
// Computing the fields
{
// Checking if we are on axis (phi is not defined there)
if (r2 == 0)
{
phi = 0;
}

cosphi = cos(phi);
cos2phi = cos(2 * phi);
sin2phi = sin(2 * phi);

// Calcultaes the E fields, assuming X polarization.
Ex_R = temp_int_i0_I + temp_int_i2_I * cos2phi;

Ex_I = -temp_int_i0_R - temp_int_i2_R * cos2phi;

Ey_R = temp_int_i2_I * sin2phi;

Ey_I = -temp_int_i2_R * sin2phi;

Ez_R = -2 * temp_int_i1_R * cosphi;

Ez_I = -2 * temp_int_i1_I * cosphi;

// Computing the amplitude

Ex = sqrt(Ex_R * Ex_R + Ex_I * Ex_I);

Ey = sqrt(Ey_R * Ey_R + Ey_I * Ey_I);

Ez = sqrt(Ez_R * Ez_R + Ez_I * Ez_I);

Ex *= PI / wavelength;
Ey *= PI / wavelength;
Ez *= PI / wavelength;

I_ = Ex_R * Ex_R
+ Ex_I
* Ex_I
+ Ey_R
* Ey_R
+ Ey_I
* Ey_I
+ Ez_R
* Ez_R
+ Ez_I
* Ez_I;

// fprintf (outFile, "%.6e\n", I_);

I_fos.writeDouble(I_);

Ex_fos.writeDouble(Ex);

Ey_fos.writeDouble(Ey);

Ez_fos.writeDouble(Ez);
}
}

private static int cacheJ(int Xsize,
int Ysize,
double scale_x,
double scale_y,
final int IntSize,
double Xc,
double Yc,
double kw,
double AlphaStart,
double theta_step,
TDoubleArrayList TheR,
TDoubleArrayList J0n,
TDoubleArrayList J1n,
TDoubleArrayList J2n,
int numOfR,
boolean found)
{
double r2;
int i;
int j;
int l;
int m;
double arg_BJ;
double temp_J0;
double temp_J1;

System.out.format("Caching Jn...\n");

for (j = 0; j < Ysize; j++)
{
for (i = 0; i < Xsize; i++)
{

// Calculate the radius
r2 = sqrt(((i - Xsize / 2 + Xc)* scale_x) * ((i - Xsize / 2 + Xc)* scale_x)
+ ((j - Ysize / 2 + Yc)* scale_y)
* ((j - Ysize / 2 + Yc)* scale_y)) ;

// Check if this radius was already calculated
for (l = 0; l < numOfR; l++)
{
if ((double) (r2 - TheR.get(l)) == 0)
{
// System.out.format("Debug: found r2=%e at l=%d, difference is %e\n",r2,
// l, r2-TheR[l]);
found = true;
break;
}

}

// if((found=false)&&(l==numOfR)){
// System.out.format("Debug: loop exit with l=%d and numOfR=%d\n", l,
// numOfR);
// }
if (!found)
{
// Radius was not found, realloc

TheR.add(r2);

for (m = 0; m < IntSize; m++)
{

arg_BJ = kw * r2 * sin(m * theta_step + AlphaStart);

// Approximate expansion
// J0n.add(BeJ(arg_BJ, 0, 50));
// J1n.add(BeJ(arg_BJ, 1, 50));
// J2n.add(BeJ(arg_BJ, 2, 50));

// Approximate polynomial expansion
temp_J0 = PolApprox_J(arg_BJ, 0);
temp_J1 = PolApprox_J(arg_BJ, 1);

J0n.add(temp_J0);
J1n.add(temp_J1);

// Calculate J2 using recurrence
if (arg_BJ == 0)
{
// Point is on the axis
J2n.add(0);
}
else
{
// Point is not on the axis
J2n.add(2 / arg_BJ * temp_J1-temp_J0);
}

}

numOfR++;

}

found = false;
l = 0;

}
}
return numOfR;
}

}