#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; }