/*
sample_zbf_program.cpp

Written by Ken Moore 5-8-2002 to demonstrate creation of a ZEMAX Beam File (ZBF)

This file is a Gaussian beam at the waist defined by parameters included below.

*/

#include "stdafx.h"
#include "stdio.h"
#include "string.h"
#include "math.h"
#include "malloc.h"
#define PI 3.14159265358979323846

int main(int argc, char* argv[])
{
int i, j, k, nx, ny, numbytes, is_polarized, unittype, xc, yc;
double x, y, dx, dy, waist, position, lambda, rayleigh, ficl, *cax, *cay, peakirr, amp, A, G, wx, wy;
FILE *out;
char filename[300];

printf("Creating sample ZBF file...\n");

/* set some parameters to use for the beam */
nx = 128;
ny = 128;
dx = 1.0; //  this is in mm if unittype is 0
dy = 1.0;
is_polarized = 0; // use 1 for a polarized beam
unittype = 0; // millimeters
waist = 12.0;
position = 0.0;
lambda = 0.6328E-03;
rayleigh = PI*waist*waist/lambda;
ficl = 0.0; //  unused for this ZBF
peakirr = 1.0;

printf("Allocating memory for beam...\n");

/* now make the beam */
numbytes = sizeof(double)*(2*(nx*ny)+1);
cax = (double *) malloc(numbytes);
cay = NULL;

if (!cax)
	{
	printf("Not enough memory for beam, aborting!\n");
	goto abort_here;
	}

if (is_polarized)
	{
	cay = (double *) malloc(numbytes);
	if (!cay)
		{
		printf("Not enough memory for beam, aborting!\n");
		goto abort_here;
		}
	}

printf("Putting Gaussian data in beam...\n");

if (is_polarized)
	{
	amp = sqrt(0.5*peakirr);
	}
else
	{
	amp = sqrt(peakirr);
	}

wx = waist;
wy = waist;
xc = 1+nx/2;
yc = 1+ny/2;
k = 1;
for (j = 1; j <= ny; j++)
	{
	for (i = 1; i <= nx; i++)
		{	
		x = (i - xc)*dx;
		y = (j - yc)*dy;
		G = (x/wx)*(x/wx) + (y/wy)*(y/wy);
		A = exp(-G);
		cax[k  ] = amp*A; // real part
		cax[k+1] = 0.0;   // imaginary part
		if (is_polarized)
			{
			cay[k  ] = amp*A;// real part
			cay[k+1] = 0.0;// imaginary part
			}
		k += 2;
		}
	}

printf("Opening beam file \"sample_beam.ZBF\" for writing...\n");

strcpy(filename,"sample_beam.ZBF");

out = fopen(filename, "wb");

// THIS IS WHERE THE BEAM DATA GETS WRITTEN OUT

/* write the format version number */
i = 0;
fwrite(&i, sizeof(int), 1, out);

/* write the beam nx size */
fwrite(&nx, sizeof(int), 1, out);

/* write the beam ny size */
fwrite(&ny, sizeof(int), 1, out);

/* write the ispolarized flag */
fwrite(&is_polarized, sizeof(int), 1, out);

/* write out the current units */
fwrite(&unittype, sizeof(int), 1, out);

/* 4 unused integers */
i = 0;
fwrite(&i, sizeof(int), 1, out);
fwrite(&i, sizeof(int), 1, out);
fwrite(&i, sizeof(int), 1, out);
fwrite(&i, sizeof(int), 1, out);

/* dx */
fwrite(&dx, sizeof(double), 1, out);

/* dy */
fwrite(&dy, sizeof(double), 1, out);

/* beam parameters */
fwrite(&position, sizeof(double), 1, out);
fwrite(&rayleigh, sizeof(double), 1, out);
fwrite(&lambda,   sizeof(double), 1, out);
fwrite(&waist,    sizeof(double), 1, out);

/* fiber coupling */
fwrite(&ficl, sizeof(double), 1, out);

/* 4 unused doubles */
x = 0.0;
fwrite(&x, sizeof(double), 1, out);
fwrite(&x, sizeof(double), 1, out);
fwrite(&x, sizeof(double), 1, out);
fwrite(&x, sizeof(double), 1, out);

/* now the beam itself */
fwrite(&cax[1], sizeof(double), nx*ny*2, out);
if (is_polarized) fwrite(&cay[1], sizeof(double), nx*ny*2, out);

fclose(out);

abort_here:;

printf("All done!\n\n\n");

if (cax) free (cax);
if (cay) free (cay);

return 0;
}




