C code for Bessel beam simulation

#include
#include
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_bessel.h>

#define intSize 512

int main (int argc, const char * argv[]) {

// argv[1] = XY size in pixel
// argv[2] = Z size in pixel
// argv[3] = wavelength in nm
// argv[4] = scale in nm/pixel
// argv[5] = NA
// argv[6] = epsilon (0 = full circular aperture, otherwise defines tha annular internal radius in terms of the NA)

if(argc != 7){
printf("Usage: %s XYsize Zsize Wavelength Scale NA epsilon\nExample: %s 100 500 488 100 0.8 0.9\n", argv[0], argv[0]);
return 1;
}

int XYsize = atoi(argv[1]);
int Zsize = atoi(argv[2]);
double wavelength = atof(argv[3]);
double scale = atof(argv[4]);

double NA = atof(argv[5]);
double epsilon = atof(argv[6]);

double Xc=0;
double Yc=0;
double Zc=0;

double RefIndex = 1.333; // Ref. index of water

double kw=2*M_PI/wavelength; // Wave vector
double MaxAlpha = asin(NA/RefIndex); // Max. half-angle defined by NA

double r2;

double z2;

double phi;

int i, j, k, l, m;

double result_array[3];
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;

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

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

// The real and imaginary parts of the electric fields
double Ex_R, Ex_I, Ey_R, Ey_I, Ez_R, Ez_I;

// The fields
double I_;
double Ex;
double Ey;
double Ez;

// Cahing
double *TheR;
double *J0n;
double *J1n;
double *J2n;

int numOfR=0;
int found;

TheR=malloc(0*sizeof(double));
if (TheR==NULL) {
printf("Error allocating memory for TheR!\n"); //print an error message
return 1; //return with failure
}

J0n=malloc(0*intSize*sizeof(double));
if (J0n==NULL) {
printf("Error allocating memory for J0n!\n"); //print an error message
return 1; //return with failure
}

J1n=malloc(0*intSize*sizeof(double));
if (J1n==NULL) {
printf("Error allocating memory for J1n!\n"); //print an error message
return 1; //return with failure
}

J2n=malloc(0*intSize*sizeof(double));
if (J2n==NULL) {
printf("Error allocating memory for J2n!\n"); //print an error message
return 1; //return with failure
}

printf("Using lambda=%d, NA=%2.2f, e=%2.2f\n", (int) wavelength, NA, epsilon);

printf("Caching Jn...\n");

for(j=0; j<XYsize; j++){
for(i=0; i<XYsize; i++){

// Calculate the radius
r2=sqrt(((i-XYsize/2+Xc))*((i-XYsize/2+Xc))+((j-XYsize/2+Yc))*((j-XYsize/2+Yc)))*scale;

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

}

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

double *temp = realloc((double*) TheR, (sizeof(&TheR))+sizeof(double) );
if ( temp != NULL ) //realloc was successful
{
TheR = temp;
TheR[numOfR] = r2;
}else{
printf("Error re-allocating memory for TheR at l=%d!\n", l); //print an error message

free(TheR);
free(J0n);
free(J1n);
free(J2n);

return 1; //return with failure
}

temp = realloc((double*) J0n, (sizeof(&J0n))+intSize*sizeof(double) );
if ( temp != NULL ) //realloc was successful
{
J0n = temp;
}else{
printf("Error re-allocating memory for J0n at l=%d!\n", l); //print an error message

free(TheR);
free(J0n);
free(J1n);
free(J2n);

return 1; //return with failure
}

temp = realloc((double*) J1n, (sizeof(&J1n))+intSize*sizeof(double) );
if ( temp != NULL ) //realloc was successful
{
J1n = temp;
}else{
printf("Error re-allocating memory for J0n at l=%d!\n", l); //print an error message

free(TheR);
free(J0n);
free(J1n);
free(J2n);

return 1; //return with failure
}

temp = realloc((double*) J2n, (sizeof(&J2n))+intSize*sizeof(double) );
if ( temp != NULL ) //realloc was successful
{
J2n = temp;
}else{
printf("Error re-allocating memory for J0n at l=%d!\n", l); //print an error message

free(TheR);
free(J0n);
free(J1n);
free(J2n);

return 1; //return with failure
}

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

// Faster (compute all three Bessel functions at once)
//gsl_sf_bessel_Jn_array(0, 2, kw*r2*sin(m*theta_step+AlphaStart), result_array);

// Slower, more accurate
result_array[0]=gsl_sf_bessel_J0(kw*r2*sin(m*theta_step+AlphaStart));
result_array[1]=gsl_sf_bessel_J1(kw*r2*sin(m*theta_step+AlphaStart));
result_array[2]=gsl_sf_bessel_Jn(2, kw*r2*sin(m*theta_step+AlphaStart));

J0n[numOfR*intSize+m] = result_array[0];

J1n[numOfR*intSize+m] = result_array[1];

J2n[numOfR*intSize+m] = result_array[2];

}

numOfR++;

}

found = 0;
l=0;

}
}

printf ("Caching finished (%d radii at %d steps)\n", numOfR, intSize);

/*
FILE *debugFile0 = fopen("Debug_r2.txt", "wt");

for(k=0; k<numOfR; k++) {

fprintf (debugFile0, "%.6e\n", TheR[k]);

}

fclose(debugFile0);

FILE *debugFile1 = fopen("Debug_J0.txt", "wt");

for(k=0; k<numOfR*intSize; k++) {

fprintf (debugFile1, "%f\n", J0n[k]);
//fwrite(&J0n[k], sizeof(double), 1, debugFile1);

}

fclose(debugFile1);
*/

char I_fileName[21];
char Ex_fileName[22];
char Ey_fileName[22];
char Ez_fileName[22];

sprintf(I_fileName, "%3d_3D_I_%2.1f_%2.1f.dat\0", (int) wavelength, NA, epsilon);

sprintf(Ex_fileName, "%3d_3D_Ex_%2.1f_%2.1f.dat\0", (int) wavelength ,NA, epsilon);

sprintf(Ey_fileName, "%3d_3D_Ey_%2.1f_%2.1f.dat\0", (int) wavelength, NA, epsilon);

sprintf(Ez_fileName, "%3d_3D_Ez_%2.1f_%2.1f.dat\0", (int) wavelength, NA, epsilon);

FILE *outFile_I;

outFile_I = fopen(I_fileName, "wb");

FILE *outFile_Ex;

outFile_Ex = fopen(Ex_fileName, "wb");

FILE *outFile_Ey;

outFile_Ey = fopen(Ey_fileName, "wb");

FILE *outFile_Ez;

outFile_Ez = fopen(Ez_fileName, "wb");

for (k=0; k<Zsize; k++) {

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

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

// Computing the position
r2=sqrt(((i-XYsize/2+Xc))*((i-XYsize/2+Xc))+((j-XYsize/2+Yc))*((j-XYsize/2+Yc)))*scale;

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

phi=asin((j-XYsize/2+Yc)*scale/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 according to r2

for(m=0; m<numOfR; m++){
if((r2-TheR[m])==0){
// Debug
//printf("Debug: difference is %f\n", (r2-TheR[m]));
break;
}

}

// Debug
//printf("Debug: r2 was %f, chosen %f\n", r2, TheR[m]);

for(l=0; l<intSize; l++){

// Computing integrals
//temp_int_common_R= theta_step*sqrt(cos(l*theta_step+AlphaStart))*cos(PhaseMask(l*theta_step+AlphaStart,phi)-kw*z2*cos(l*theta_step+AlphaStart))*sin(l*theta_step+AlphaStart);
//temp_int_common_I=theta_step*sqrt(cos(l*theta_step+AlphaStart))*sin(PhaseMask(l*theta_step+AlphaStart,phi)-kw*z2*cos(l*theta_step+AlphaStart))*sin(l*theta_step+AlphaStart);
temp_int_common_R= theta_step*sqrt(cos(l*theta_step+AlphaStart))*cos(kw*z2*cos(l*theta_step+AlphaStart))*sin(l*theta_step+AlphaStart);
temp_int_common_I=theta_step*sqrt(cos(l*theta_step+AlphaStart))*sin(-kw*z2*cos(l*theta_step+AlphaStart))*sin(l*theta_step+AlphaStart);

temp_int_i0_R+=temp_int_common_R*J0n[m*intSize+l]*(1+cos(l*theta_step+AlphaStart));

temp_int_i0_I+=temp_int_common_I*J0n[m*intSize+l]*(1+cos(l*theta_step+AlphaStart));

temp_int_i1_R+=temp_int_common_R*J1n[m*intSize+l]*sin(l*theta_step+AlphaStart);

temp_int_i1_I+=temp_int_common_I*J1n[m*intSize+l]*sin(l*theta_step+AlphaStart);

temp_int_i2_R+=temp_int_common_R*J2n[m*intSize+l]*(1-cos(l*theta_step+AlphaStart));

temp_int_i2_I+=temp_int_common_I*J2n[m*intSize+l]*(1-cos(l*theta_step+AlphaStart));

}

// Computing the fields

// Checking if we are on axis (phi is not defined there)
if(r2==0){
phi = 0;
}

Ex_R=-temp_int_i0_R+temp_int_i2_R*cos(2*phi);

Ex_I=temp_int_i0_I+temp_int_i2_I*cos(2*phi);

Ey_R=-temp_int_i2_R*sin(2*phi);

Ey_I=temp_int_i2_I*sin(2*phi);

Ez_R=-2*temp_int_i1_R*cos(phi);

Ez_I=-2*temp_int_i1_I*cos(phi);

// Computing the amplitude

Ex=sqrt(Ex_R*Ex_R+Ex_I*Ex_I)*M_PI/wavelength;

Ey=sqrt(Ey_R*Ey_R+Ey_I*Ey_I)*M_PI/wavelength;

Ez=sqrt(Ez_R*Ez_R+Ez_I*Ez_I)*M_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_);

fwrite(&I_, sizeof(double), 1, outFile_I);

fwrite(&Ex, sizeof(double), 1, outFile_Ex);

fwrite(&Ey, sizeof(double), 1, outFile_Ey);

fwrite(&Ez, sizeof(double), 1, outFile_Ez);

}

}

//printf("Computing plane %d of %d\n", k+1, Zsize);
}

fclose (outFile_I);
fclose (outFile_Ex);
fclose (outFile_Ey);
fclose (outFile_Ez);

free(TheR);
free(J0n);
free(J1n);
free(J2n);

return 0;
}