Diffraction.c
/*
http://www.stfmc.de/misc/diffcontrarefl/Diffraction.c
belonging to
http://www.stfmc.de/misc/diffcontrarefl/tlf.html
Light sources with spikes.
Apply a rigorous diffraction theory to a camera situation.
The numerical part.
This program was written by Stephan Seidl in 2014, no copyright is claimed.
It is offered as-is, without any warranty.
This program is in the public domain; do with it what you wish.
This is version 1.1 of the program.
The names of the variables recvrrank and sendrrank
have been exchanged for a better understanding.
This does not affect the code.
*/
#ifndef _GNU_SOURCE
#define _GNU_SOURCE
#endif
#include <stdio.h>
#include <sys/time.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#if (0)
#ifndef MYDEBUG
#define MYDEBUG
#endif
#endif
#if (defined (MYRUNATALMPS1) || defined (MYRUNATALMPS0))
/* Asymptotically Large Massively Parallel System */
#ifndef MYUSEMPI
#define MYUSEMPI
#endif
#endif
#ifdef MYUSEMPI
#include <mpi.h>
#ifdef MYDEBUG
#undef MYDEBUG
#endif
#endif
static void heapsortsrcpoints (double *, unsigned int *, unsigned int *, unsigned int);
static double x2nf (const char *);
int main (int argc, char **argv)
{
double waveLengthTable [10] [5];
char buf [1000];
int isize, irank;
unsigned int it, jt, idiaphragmpos, isrcdisplacement, kwavelength,
nsrcpoints, nblades, nrad, irowmin, irowbeyond, jcolmin, jcolbeyond,
uselargetables, iwavelengthmin, iwavelengthbeyond, ip, jp, *rankbypixno,
npixperrank, rankentry, ipixperrank, *printbypixno, npixperprint,
wkunitentry, ipixperprint, *srcpointindexrow, *srcpointindexcol,
nsrcpointsmax, is, js, *nphibyirad, *tableindexccsbyirad, sumnphi, irad,
nphi, nphirem, megstomalloc, indexccsbyirad, iphi, indexccs, megsmalloced,
isrcpoint, iwavelength, numrankbits, kauxrankbits, ibit, commbitmask,
garbagemask, recvrrank, sendrrank, inormalization, pixr, pixg, pixb;
double focalLengthObjective, fNumber, distanceSourceLensLeft,
distanceLensLeftLensRight, streetLightImageOffsetHorizontal,
streetLightImageOffsetVertical, streetLightImageMajorSemiAxis,
streetLightImageMinorSemiAxis, twoxpi, fourxpireci, diameterEntrancePupil,
distanceLensLeftDiaphragm, distanceDiaphragmLensRight,
focalLengthLensLeft, focalLengthLensRight, distanceLensRightTarget,
phiblade, auxsin, auxcos, auxnum, auxden,
radiusCircumCircleForDiameterTwoEquiv, radmaxmax, drad, rsz, rdz, rtz,
rssz, rttz, rsssz, rtttz, rssdz, rttdz, *picturer, *pictureg, *pictureb,
*srcpointdistance, a, b, e, rsyseeming, rsxseeming, aux2a, rad,
*tableccsphi, dphi, phi, phirem, sinphi, cosphi, wkunitsecs, rtxseeming,
rtyseeming, rtx, rty, rttx, rtty, rtttx, rttty, rtttabs, rtttcos, rsx,
rsy, rssx, rssy, rsssx, rsssy, rsssabs, rssscos, k, pixelaccureal,
pixelaccuimag, radmax, differential, rdx, rdy, rssdx, rssdy, deltas,
rttdx, rttdy, deltat, paths, patht, gampli0, gampli1, gampli2, angle,
sinangle, cosangle, intensity, *picturealt, pixc1old, pixc0old, vmax,
vmin, picscale, pixstandarddeviation, pixaverage, pixdeviation, pixc1,
pixc0, pixvalmapped, epsiter, apixr, apixg, apixb;
struct timeval tv0, tv1;
#if (defined (MYUSEMPI))
MPI_Status status;
#endif
FILE *pfile;
#if (defined (MYUSEMPI))
(void) MPI_Init (&argc, &argv);
(void) MPI_Comm_size ((MPI_Comm) MPI_COMM_WORLD, &isize);
(void) MPI_Comm_rank ((MPI_Comm) MPI_COMM_WORLD, &irank);
#else
isize = 1;
irank = 0;
#endif
/* Strong lines taken from http://physics.nist.gov/PhysRefData/Handbook/element_name.htm;
RGB values taken from http://www.magnetkern.de/spektrum.html,
wavelength column in accordance with CIECAM02. */
waveLengthTable[0][0] = 643.8e-9; /* Cadmium Cd I */
waveLengthTable[0][1] = x2nf ("ff");
waveLengthTable[0][2] = x2nf ("00");
waveLengthTable[0][3] = x2nf ("1c");
waveLengthTable[1][0] = 589.3e-9; /* Sodium D1/D2 doublet */
waveLengthTable[1][1] = x2nf ("ff");
waveLengthTable[1][2] = x2nf ("a0");
waveLengthTable[1][3] = x2nf ("00");
waveLengthTable[2][0] = 546.1e-9; /* Mercury Hg I */
waveLengthTable[2][1] = x2nf ("56");
waveLengthTable[2][2] = x2nf ("ff");
waveLengthTable[2][3] = x2nf ("00");
waveLengthTable[3][0] = 484.4e-9; /* Xenon Xe II */
waveLengthTable[3][1] = x2nf ("00");
waveLengthTable[3][2] = x2nf ("e5");
waveLengthTable[3][3] = x2nf ("ff");
waveLengthTable[4][0] = 451.1e-9; /* Indium In I */
waveLengthTable[4][1] = x2nf ("2e");
waveLengthTable[4][2] = x2nf ("00");
waveLengthTable[4][3] = x2nf ("ff");
waveLengthTable[0][4] = 0.6304558129449025;
waveLengthTable[1][4] = 0.6638371029589381;
waveLengthTable[2][4] = 0.6950309656576327;
waveLengthTable[3][4] = 0.8794840815443867;
waveLengthTable[4][4] = 1.0000000000000000;
for (it = 0; it < 5; it++) {
waveLengthTable[5 + it][0] = 20.0 * waveLengthTable[it][0];
for (jt = 1; jt < 5; jt++)
waveLengthTable[5 + it][jt] = waveLengthTable[it][jt];
}
focalLengthObjective = 23.0e-3;
fNumber = 8.0;
distanceSourceLensLeft = 80.0;
distanceLensLeftLensRight = 30.0e-3;
streetLightImageOffsetHorizontal = 1.53e-3;
streetLightImageOffsetVertical = 3.67e-3;
streetLightImageMajorSemiAxis = 0.05e-3;
streetLightImageMinorSemiAxis = 0.03e-3;
#if (defined (MYDEBUG))
idiaphragmpos = 2;
isrcdisplacement = 5;
kwavelength = 1;
nsrcpoints = 1;
nblades = 3;
nrad = 100;
irowmin = 0;
irowbeyond = 400;
jcolmin = 0;
jcolbeyond = 600;
uselargetables = 1;
#elif (defined(MYRUNATALMPS1))
/* Assume 160 hours with 240000 P4@2.5GHz processors. */
idiaphragmpos = 2;
isrcdisplacement = 5;
kwavelength = 10;
nsrcpoints = 157;
nblades = 7;
nrad = 20000;
irowmin = 0;
irowbeyond = 400;
jcolmin = 0;
jcolbeyond = 600;
uselargetables = 0;
#elif (defined(MYRUNATALMPS0))
/* Assume 1 hour with 240000 P4@2.5GHz processors. */
idiaphragmpos = 2;
isrcdisplacement = 5;
kwavelength = 10;
nsrcpoints = 1;
nblades = 7;
nrad = 20000;
irowmin = 0;
irowbeyond = 400;
jcolmin = 0;
jcolbeyond = 600;
uselargetables = 0;
#else
(void) strcpy (&buf[0], "Usage: ");
(void) strcat (&buf[0], *(argv + 0));
(void) strcat (&buf[0], " <idiaphragmpos>");
(void) strcat (&buf[0], " <isrcdisplacement>");
(void) strcat (&buf[0], " <kwavelength>");
(void) strcat (&buf[0], " \\\n ");
(void) strcat (&buf[0], " <nsrcpoints>");
(void) strcat (&buf[0], " <nblades>");
(void) strcat (&buf[0], " <nrad>");
(void) strcat (&buf[0], " <irowmin>");
(void) strcat (&buf[0], " <irowbeyond>");
(void) strcat (&buf[0], " \\\n ");
(void) strcat (&buf[0], " <jcolmin>");
(void) strcat (&buf[0], " <jcolbeyond>");
(void) strcat (&buf[0], " <uselargetables>");
(void) strcat (&buf[0], "\n");
if (argc != 12) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 1), "%u", &idiaphragmpos) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 2), "%u", &isrcdisplacement) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 3), "%u", &kwavelength) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 4), "%u", &nsrcpoints) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 5), "%u", &nblades) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 6), "%u", &nrad) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 7), "%u", &irowmin) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 8), "%u", &irowbeyond) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 9), "%u", &jcolmin) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 10), "%u", &jcolbeyond) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
if (sscanf (*(argv + 11), "%u", &uselargetables) != 1) {
(void) printf ("%s", &buf[0]);
return (1);
}
#endif
if (idiaphragmpos > 8)
idiaphragmpos = 8;
if (isrcdisplacement > 5)
isrcdisplacement = 5;
if (kwavelength > 11)
kwavelength = 11;
if (nsrcpoints < 1)
nsrcpoints = 1;
if (nsrcpoints > 999)
nsrcpoints = 999;
if (nblades < 3)
nblades = 3;
if (nblades > 19)
nblades = 19;
if (nrad < 1)
nrad = 1;
if (nrad > 99999)
nrad = 99999;
if (irowmin > 399)
irowmin = 399;
if (irowbeyond < irowmin + 1)
irowbeyond = irowmin + 1;
if (irowbeyond > 400)
irowbeyond = 400;
if (jcolmin > 599)
jcolmin = 599;
if (jcolbeyond < jcolmin + 1)
jcolbeyond = jcolmin + 1;
if (jcolbeyond > 600)
jcolbeyond = 600;
if (uselargetables > 1)
uselargetables = 1;
if (kwavelength == 10) {
iwavelengthmin = 0;
iwavelengthbeyond = 5;
}
else if (kwavelength == 11) {
iwavelengthmin = 5;
iwavelengthbeyond = 10;
}
else {
iwavelengthmin = kwavelength;
iwavelengthbeyond = kwavelength + 1;
}
twoxpi = 8.0 * atan (1.0);
fourxpireci = 1.0 / (2.0 * twoxpi);
diameterEntrancePupil = focalLengthObjective / fNumber;
distanceLensLeftDiaphragm = distanceLensLeftLensRight * ((double) idiaphragmpos * 0.125);
distanceDiaphragmLensRight = distanceLensLeftLensRight - distanceLensLeftDiaphragm;
focalLengthLensLeft = distanceSourceLensLeft;
focalLengthLensRight = focalLengthObjective * (focalLengthLensLeft - distanceLensLeftLensRight)
/ (focalLengthLensLeft - focalLengthObjective);
distanceLensRightTarget = focalLengthLensRight;
phiblade = twoxpi / (double) nblades;
(void) sincos (phiblade * 0.5, &auxsin, &auxcos);
auxnum = sqrt ((double) nblades * auxsin * auxcos * auxcos * auxcos);
auxden = (double) nblades * auxsin * auxcos * auxcos;
radiusCircumCircleForDiameterTwoEquiv = sqrt (twoxpi * 0.5) * auxnum / auxden;
radmaxmax = radiusCircumCircleForDiameterTwoEquiv * diameterEntrancePupil / 2.0;
drad = radmaxmax / (double) nrad;
rsz = -(distanceSourceLensLeft + distanceLensLeftDiaphragm);
rdz = 0.0;
rtz = distanceDiaphragmLensRight + distanceLensRightTarget;
rssz = rdz;
rttz = rdz;
rsssz = rssz - rsz;
rtttz = rtz - rttz;
rssdz = rssz - rdz;
rttdz = rttz - rdz;
picturer = (double *) malloc (sizeof (double) * 400 * 600);
if (picturer == (double *) 0) {
(void) printf ("No more core.\n");
return (1);
}
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(picturer + ip * 600 + jp) = 0.0;
pictureg = (double *) malloc (sizeof (double) * 400 * 600);
if (pictureg == (double *) 0) {
(void) printf ("No more core.\n");
return (1);
}
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(pictureg + ip * 600 + jp) = 0.0;
pictureb = (double *) malloc (sizeof (double) * 400 * 600);
if (pictureb == (double *) 0) {
(void) printf ("No more core.\n");
return (1);
}
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(pictureb + ip * 600 + jp) = 0.0;
rankbypixno = (unsigned int *) malloc (sizeof (unsigned int) * 400 * 600);
if (rankbypixno == (unsigned int *) 0) {
(void) printf ("No more core.\n");
return (1);
}
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(rankbypixno + ip * 600 + jp) = (unsigned int) ~ (unsigned int) 0;
npixperrank = (irowbeyond - irowmin) * (jcolbeyond - jcolmin) / (unsigned int) isize;
if (npixperrank * (unsigned int) isize != (irowbeyond - irowmin) * (jcolbeyond - jcolmin))
npixperrank++;
rankentry = 0;
ipixperrank = 0;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
*(rankbypixno + ip * 600 + jp) = rankentry;
ipixperrank++;
if (ipixperrank >= npixperrank) {
rankentry++;
ipixperrank = 0;
}
}
printbypixno = (unsigned int *) malloc (sizeof (unsigned int) * 400 * 600);
if (printbypixno == (unsigned int *) 0) {
(void) printf ("No more core.\n");
return (1);
}
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(printbypixno + ip * 600 + jp) = 0;
npixperprint = npixperrank / 1000;
if (npixperprint * 1000 != npixperrank)
npixperprint++;
wkunitentry = 0;
ipixperprint = 0;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
if (*(rankbypixno + ip * 600 + jp) != 0)
break;
if (ip != irowmin || jp != jcolmin)
ipixperprint++;
if (ipixperprint >= npixperprint) {
*(printbypixno + ip * 600 + jp) = ++wkunitentry;
ipixperprint = 0;
}
}
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
if (*(printbypixno + ip * 600 + jp) == 0)
continue;
*(printbypixno + ip * 600 + jp) = wkunitentry - *(printbypixno + ip * 600 + jp);
}
srcpointdistance = (double *) malloc (sizeof (double) * 400 * 600);
if (srcpointdistance == (double *) 0) {
(void) printf ("No more core.\n");
return (1);
}
srcpointindexrow = (unsigned int *) malloc (sizeof (unsigned int) * 400 * 600);
if (srcpointindexrow == (unsigned int *) 0) {
(void) printf ("No more core.\n");
return (1);
}
srcpointindexcol = (unsigned int *) malloc (sizeof (unsigned int) * 400 * 600);
if (srcpointindexcol == (unsigned int *) 0) {
(void) printf ("No more core.\n");
return (1);
}
nsrcpointsmax = 0;
a = streetLightImageMajorSemiAxis;
b = streetLightImageMinorSemiAxis;
e = sqrt (a * a - b * b);
for (is = 0; is < 400; is++) {
rsyseeming = (((double) is - 200.0) / 200.0) * 1.1e-3;
for (js = 0; js < 600; js++) {
rsxseeming = (((double) js - 300.0) / 300.0) * 1.65e-3;
aux2a = sqrt (rsyseeming * rsyseeming + (rsxseeming + e) * (rsxseeming + e))
+ sqrt (rsyseeming * rsyseeming + (rsxseeming - e) * (rsxseeming - e));
if (aux2a > 2.0 * a)
continue;
*(srcpointdistance + nsrcpointsmax) = sqrt (rsxseeming * rsxseeming + rsyseeming * rsyseeming);
*(srcpointindexrow + nsrcpointsmax) = is;
*(srcpointindexcol + nsrcpointsmax) = js;
nsrcpointsmax++;
}
}
if (nsrcpointsmax == 0) {
*(srcpointdistance + nsrcpointsmax) = 0.0;
*(srcpointindexrow + nsrcpointsmax) = 200;
*(srcpointindexcol + nsrcpointsmax) = 300;
nsrcpointsmax++;
}
(void) heapsortsrcpoints (srcpointdistance, srcpointindexrow, srcpointindexcol, nsrcpointsmax);
if (nsrcpoints > nsrcpointsmax)
nsrcpoints = nsrcpointsmax;
nphibyirad = (unsigned int *) malloc (sizeof (unsigned int) * (size_t) nrad);
if (nphibyirad == (unsigned int *) 0) {
(void) printf ("No more core.\n");
return (1);
}
if (uselargetables != 0) {
tableindexccsbyirad = (unsigned int *) malloc (sizeof (unsigned int) * (size_t) nrad);
if (tableindexccsbyirad == (unsigned int *) 0) {
(void) printf ("No more core.\n");
return (1);
}
}
else
tableindexccsbyirad = (unsigned int *) 0;
sumnphi = 0;
for (irad = 0; irad < nrad; irad++) {
rad = drad * ((double) irad + 0.5);
nphi = (unsigned int) (ceil (twoxpi * rad / drad) + 0.5);
if (nphi == 0) {
(void) printf ("Unrecoverable error.\n");
return (1);
}
nphirem = nphi % (2 * nblades);
if (nphirem != 0)
nphi += 2 * nblades - nphirem;
*(nphibyirad + irad) = nphi;
if (uselargetables == 0)
continue;
if (sumnphi == 0) {
*(tableindexccsbyirad + irad) = 0;
sumnphi = nphi;
continue;
}
if (*(nphibyirad + irad - 0) == *(nphibyirad + irad - 1)) {
*(tableindexccsbyirad + irad - 0) = *(tableindexccsbyirad + irad - 1);
continue;
}
*(tableindexccsbyirad + irad) = sumnphi;
sumnphi += nphi;
}
if (irank == 0) {
(void) printf ("Running: %s %u %u %u %u %u %u %u %u %u %u %u\n", *(argv + 0),
idiaphragmpos, isrcdisplacement, kwavelength, nsrcpoints, nblades,
nrad, irowmin, irowbeyond, jcolmin, jcolbeyond, uselargetables);
if (uselargetables != 0)
megstomalloc = (unsigned int) (( (double) (400 * 600 * 4 + sumnphi * 3) *
(double) sizeof (double)
+ (double) (400 * 600 * 4 + nrad * 2) *
(double) sizeof (unsigned int)) / 1.048576e+6 + 0.5);
else
megstomalloc = (unsigned int) (( (double) (400 * 600 * 4) *
(double) sizeof (double)
+ (double) (400 * 600 * 4 + nrad) *
(double) sizeof (unsigned int)) / 1.048576e+6 + 0.5);
if (megstomalloc == 0)
megstomalloc = 1;
(void) printf ("Info: %s is going to malloc %u MiB per process\n", *(argv + 0), megstomalloc);
}
if (uselargetables != 0) {
tableccsphi = (double *) malloc (sizeof (double) * (size_t) sumnphi * 3);
if (tableccsphi == (double *) 0) {
(void) printf ("No more core.\n");
return (1);
}
for (irad = 0; irad < nrad; irad++) {
rad = drad * ((double) irad + 0.5);
nphi = *(nphibyirad + irad);
dphi = twoxpi / (double) nphi;
indexccsbyirad = *(tableindexccsbyirad + irad);
for (iphi = 0; iphi < nphi; iphi++) {
phi = dphi * ((double) iphi + 0.5);
phirem = phi - round (phi / phiblade) * phiblade;
while (phirem < 0.0)
phirem += phiblade;
while (phirem > phiblade)
phirem -= phiblade;
if (phirem < 0.0)
phirem += phiblade;
if (phirem > phiblade * 0.5)
phirem = phiblade - phirem;
(void) sincos (phi, &sinphi, &cosphi);
indexccs = 3 * (indexccsbyirad + iphi);
*(tableccsphi + indexccs + 0) = cos (phirem);
*(tableccsphi + indexccs + 1) = cosphi;
*(tableccsphi + indexccs + 2) = sinphi;
}
}
}
else
tableccsphi = (double *) 0;
if (irank == 0) {
if (uselargetables != 0)
megsmalloced = (unsigned int) (( (double) (400 * 600 * 4 + sumnphi * 3) *
(double) sizeof (double)
+ (double) (400 * 600 * 4 + nrad * 2) *
(double) sizeof (unsigned int)) / 1.048576e+6 + 0.5);
else
megsmalloced = (unsigned int) (( (double) (400 * 600 * 4) *
(double) sizeof (double)
+ (double) (400 * 600 * 4 + nrad) *
(double) sizeof (unsigned int)) / 1.048576e+6 + 0.5);
if (megsmalloced == 0)
megsmalloced = 1;
(void) printf ("Info: %s malloced and attached %u MiB per process\n", *(argv + 0), megsmalloced);
}
tv0.tv_sec = tv1.tv_sec = 0;
tv0.tv_usec = tv1.tv_usec = 0;
for (it = irowmin; it < irowbeyond; it++)
for (jt = jcolmin; jt < jcolbeyond; jt++) {
#if (defined (MYDEBUG))
if (it != 80 || jt != 150)
#else
if (*(rankbypixno + it * 600 + jt) != (unsigned int) irank)
#endif
continue;
if (irank == 0 && it == irowmin && jt == jcolmin) {
if (gettimeofday (&tv1, (struct timezone *) 0) != 0) {
(void) printf ("gettimeofday() failed.\n");
return (1);
}
tv0.tv_sec = tv1.tv_sec;
tv0.tv_usec = tv1.tv_usec;
}
else if (*(printbypixno + it * 600 + jt) != 0) {
if (gettimeofday (&tv1, (struct timezone *) 0) != 0) {
(void) printf ("gettimeofday() failed.\n");
return (1);
}
wkunitsecs = ((double) tv1.tv_sec - (double) tv0.tv_sec )
+ ((double) tv1.tv_usec - (double) tv0.tv_usec) * 1.0e-6;
tv0.tv_sec = tv1.tv_sec;
tv0.tv_usec = tv1.tv_usec;
(void) printf ("Running %s: wkunitsecs=%.3f wkunitstodo=%u\n",
*(argv + 0), wkunitsecs, *(printbypixno + it * 600 + jt));
}
rtxseeming = (((double) jt - 300.0) / 300.0) * 1.65e-3;
rtyseeming = (((double) it - 200.0) / 200.0) * 1.1e-3;
rtx = (rtxseeming + ((double) isrcdisplacement * 0.2) * streetLightImageOffsetHorizontal);
rty = (rtyseeming + ((double) isrcdisplacement * 0.2) * streetLightImageOffsetVertical);
rttx = rtx * (-(distanceDiaphragmLensRight / distanceLensRightTarget));
rtty = rty * (-(distanceDiaphragmLensRight / distanceLensRightTarget));
rtttx = rtx - rttx;
rttty = rty - rtty;
rtttabs = sqrt (rtttx * rtttx + rttty * rttty + rtttz * rtttz);
rtttcos = rtttz / rtttabs;
for (isrcpoint = 0; isrcpoint < nsrcpoints; isrcpoint++) {
is = *(srcpointindexrow + isrcpoint);
js = *(srcpointindexcol + isrcpoint);
rsxseeming = (((double) js - 300.0) / 300.0) * 1.65e-3;
rsyseeming = (((double) is - 200.0) / 200.0) * 1.1e-3;
rsx = (rsxseeming + ((double) isrcdisplacement * 0.2) * streetLightImageOffsetHorizontal)
* (-(distanceSourceLensLeft / distanceLensRightTarget));
rsy = (rsyseeming + ((double) isrcdisplacement * 0.2) * streetLightImageOffsetVertical)
* (-(distanceSourceLensLeft / distanceLensRightTarget));
rssx = rsx * (-(distanceLensLeftDiaphragm / distanceSourceLensLeft));
rssy = rsy * (-(distanceLensLeftDiaphragm / distanceSourceLensLeft));
rsssx = rssx - rsx;
rsssy = rssy - rsy;
rsssabs = sqrt (rsssx * rsssx + rsssy * rsssy + rsssz * rsssz);
rssscos = rsssz / rsssabs;
for (iwavelength = iwavelengthmin; iwavelength < iwavelengthbeyond; iwavelength++) {
k = twoxpi / waveLengthTable[iwavelength][0];
pixelaccureal = pixelaccuimag = 0.0;
for (irad = 0; irad < nrad; irad++) {
#if (defined (MYDEBUG))
if (irad != 23)
continue;
#endif
rad = drad * ((double) irad + 0.5);
nphi = *(nphibyirad + irad);
dphi = twoxpi / (double) nphi;
if (uselargetables != 0)
indexccsbyirad = *(tableindexccsbyirad + irad);
else
indexccsbyirad = 0;
for (iphi = 0; iphi < nphi; iphi++) {
#if (defined (MYDEBUG))
if (iphi != 43)
continue;
#endif
if (uselargetables != 0) {
phi = 0.0;
indexccs = 3 * (indexccsbyirad + iphi);
radmax = radmaxmax * *(tableccsphi + indexccs + 0);
}
else {
indexccs = 0;
phi = dphi * ((double) iphi + 0.5);
phirem = phi - round (phi / phiblade) * phiblade;
while (phirem < 0.0)
phirem += phiblade;
while (phirem > phiblade)
phirem -= phiblade;
if (phirem < 0.0)
phirem += phiblade;
if (phirem > phiblade * 0.5)
phirem = phiblade - phirem;
radmax = radmaxmax * cos (phirem);
}
#if (defined (MYDEBUG))
differential = dphi * rad * drad;
#else
if (rad - 0.5 * drad > radmax)
continue;
if (rad + 0.5 * drad < radmax)
differential = dphi * rad * drad;
else
differential = dphi * rad * (radmax - (rad - 0.5 * drad));
#endif
if (uselargetables != 0) {
rdx = rad * *(tableccsphi + indexccs + 1);
rdy = rad * *(tableccsphi + indexccs + 2);
}
else {
(void) sincos (phi, &sinphi, &cosphi);
rdx = rad * cosphi;
rdy = rad * sinphi;
}
rssdx = rssx - rdx;
rssdy = rssy - rdy;
deltas = (rsssx * rssdx + rsssy * rssdy + rsssz * rssdz) / rsssabs;
rttdx = rttx - rdx;
rttdy = rtty - rdy;
deltat = (rtttx * rttdx + rttty * rttdy + rtttz * rttdz) / rtttabs;
paths = rsssabs - deltas;
patht = rtttabs + deltat;
gampli0 = paths * patht;
gampli0 = k * gampli0 * gampli0;
gampli0 = waveLengthTable[iwavelength][4] * fourxpireci * differential / gampli0;
gampli1 = gampli0 * k * paths * patht * (rssscos + rtttcos);
gampli2 = gampli0 * (paths * rtttcos + patht * rssscos);
angle = k * (paths + patht);
angle -= round (angle / twoxpi) * twoxpi;
#if (defined (MYDEBUG))
while (angle < 0.0)
angle += twoxpi;
while (angle > twoxpi)
angle -= twoxpi;
if (angle < 0.0)
angle += twoxpi;
#endif
(void) sincos (angle, &sinangle, &cosangle);
pixelaccureal += gampli1 * sinangle + gampli2 * cosangle;
pixelaccuimag += gampli1 * cosangle - gampli2 * sinangle;
#if (defined (MYDEBUG))
printf("\n");
printf("# nphi := %d:\n", (int) nphi);
printf("# k := %22.15e:\n", k);
printf("# rdx := %22.15e: rdy := %22.15e: rdz := %22.15e:\n", rdx, rdy, rdz);
printf("# rsx := %22.15e: rsy := %22.15e: rsz := %22.15e:\n", rsx, rsy, rsz);
printf("# rssx := %22.15e: rssy := %22.15e: rssz := %22.15e:\n", rssx, rssy, rssz);
printf("# rsssx := %22.15e: rsssy := %22.15e: rsssz := %22.15e:\n", rsssx, rsssy, rsssz);
printf("# rsssabs := %22.15e: rssscos := %22.15e:\n", rsssabs, rssscos);
printf("# rssdx := %22.15e: rssdy := %22.15e: rssdz := %22.15e:\n", rssdx, rssdy, rssdz);
printf("# deltasx := %22.15e: deltasy := %22.15e: deltasz := %22.15e:\n", rsssx * rssdx / rsssabs,
rsssy * rssdy / rsssabs,
rsssz * rssdz / rsssabs);
printf("# deltas := %22.15e:\n", deltas);
printf("# rtx := %22.15e: rty := %22.15e: rtz := %22.15e:\n", rtx, rty, rtz);
printf("# rttx := %22.15e: rtty := %22.15e: rttz := %22.15e:\n", rttx, rtty, rttz);
printf("# rtttx := %22.15e: rttty := %22.15e: rtttz := %22.15e:\n", rtttx, rttty, rtttz);
printf("# rtttabs := %22.15e: rtttcos := %22.15e:\n", rtttabs, rtttcos);
printf("# rttdx := %22.15e: rttdy := %22.15e: rttdz := %22.15e:\n", rttdx, rttdy, rttdz);
printf("# deltatx := %22.15e: deltaty := %22.15e: deltatz := %22.15e:\n", rtttx * rttdx / rtttabs,
rttty * rttdy / rtttabs,
rtttz * rttdz / rtttabs);
printf("# deltat := %22.15e:\n", deltat);
printf("# paths := %22.15e: patht := %22.15e: pathtot := %22.15e:\n", paths, patht, paths + patht);
printf("# gampli0 := %22.15e: gampli1 := %22.15e: gampli2 := %22.15e:\n", gampli0, gampli1, gampli2);
printf("# angle := %22.15e: cosangle := %22.15e: sinangle := %22.15e:\n", angle, cosangle, sinangle);
printf("\n");
#endif
}
}
intensity = pixelaccureal * pixelaccureal + pixelaccuimag * pixelaccuimag;
*(picturer + it * 600 + jt) += waveLengthTable[iwavelength][1] * intensity;
*(pictureg + it * 600 + jt) += waveLengthTable[iwavelength][2] * intensity;
*(pictureb + it * 600 + jt) += waveLengthTable[iwavelength][3] * intensity;
}
}
}
if (uselargetables != 0) {
(void) free (tableccsphi);
(void) free (tableindexccsbyirad);
}
(void) free (nphibyirad);
(void) free (srcpointindexcol);
(void) free (srcpointindexrow);
(void) free (srcpointdistance);
(void) free (printbypixno);
(void) free (rankbypixno);
numrankbits = 0;
kauxrankbits = (unsigned int) isize - 1;
while (kauxrankbits != 0) {
kauxrankbits >>= 1;
numrankbits++;
}
for (ibit = 0; ibit < numrankbits; ibit++) {
commbitmask = (unsigned int) 0x1 << ibit;
garbagemask = commbitmask - 1;
if (((unsigned int) irank & garbagemask) != 0)
continue;
recvrrank = (unsigned int) irank & ~ commbitmask;
sendrrank = (unsigned int) irank | commbitmask;
if (sendrrank >= (unsigned int) isize)
continue;
if ((unsigned int) irank == recvrrank) {
picturealt = (double *) malloc (sizeof (double) * 400 * 600);
if (picturealt == (double *) 0) {
(void) printf ("No more core.\n");
return (1);
}
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(picturealt + ip * 600 + jp) = 0.0;
#if (defined (MYUSEMPI))
(void) MPI_Recv ((void *) picturealt, 400 * 600 * sizeof (double), MPI_BYTE,
(int) sendrrank, 0, (MPI_Comm) MPI_COMM_WORLD, &status);
buf[0] = '\0';
(void) MPI_Send ((void *) &buf[0], sizeof (char), MPI_BYTE,
(int) sendrrank, 0, (MPI_Comm) MPI_COMM_WORLD);
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(picturer + ip * 600 + jp) += *(picturealt + ip * 600 + jp);
(void) MPI_Recv ((void *) picturealt, 400 * 600 * sizeof (double), MPI_BYTE,
(int) sendrrank, 0, (MPI_Comm) MPI_COMM_WORLD, &status);
buf[0] = '\0';
(void) MPI_Send ((void *) &buf[0], sizeof (char), MPI_BYTE,
(int) sendrrank, 0, (MPI_Comm) MPI_COMM_WORLD);
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(pictureg + ip * 600 + jp) += *(picturealt + ip * 600 + jp);
(void) MPI_Recv ((void *) picturealt, 400 * 600 * sizeof (double), MPI_BYTE,
(int) sendrrank, 0, (MPI_Comm) MPI_COMM_WORLD, &status);
buf[0] = '\0';
(void) MPI_Send ((void *) &buf[0], sizeof (char), MPI_BYTE,
(int) sendrrank, 0, (MPI_Comm) MPI_COMM_WORLD);
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
*(pictureb + ip * 600 + jp) += *(picturealt + ip * 600 + jp);
#else
;
#endif
(void) free (picturealt);
}
else if ((unsigned int) irank == sendrrank) {
#if (defined (MYUSEMPI))
(void) MPI_Send ((void *) picturer, 400 * 600 * sizeof (double), MPI_BYTE,
(int) recvrrank, 0, (MPI_Comm) MPI_COMM_WORLD);
(void) MPI_Recv ((void *) &buf[0], sizeof (char), MPI_BYTE,
(int) recvrrank, 0, (MPI_Comm) MPI_COMM_WORLD, &status);
(void) MPI_Send ((void *) pictureg, 400 * 600 * sizeof (double), MPI_BYTE,
(int) recvrrank, 0, (MPI_Comm) MPI_COMM_WORLD);
(void) MPI_Recv ((void *) &buf[0], sizeof (char), MPI_BYTE,
(int) recvrrank, 0, (MPI_Comm) MPI_COMM_WORLD, &status);
(void) MPI_Send ((void *) pictureb, 400 * 600 * sizeof (double), MPI_BYTE,
(int) recvrrank, 0, (MPI_Comm) MPI_COMM_WORLD);
(void) MPI_Recv ((void *) &buf[0], sizeof (char), MPI_BYTE,
(int) recvrrank, 0, (MPI_Comm) MPI_COMM_WORLD, &status);
#else
;
#endif
}
else {
(void) printf ("Unrecoverable error.\n");
return (1);
}
}
if (irank == 0) {
(void) sprintf (&buf[0], "DiffractionI%uD%uW%02uS%03uB%02uR%05uV%03uL%03uH%03uL%03u.data",
idiaphragmpos, isrcdisplacement, kwavelength, nsrcpoints, nblades,
nrad, irowmin, irowbeyond, jcolmin, jcolbeyond);
pfile = fopen (&buf[0], "w");
if (pfile == (FILE *) 0) {
(void) printf ("Unable to open %s.\n", &buf[0]);
return (1);
}
(void) fprintf (pfile, "P3\n");
(void) fprintf (pfile, "600 400\n");
(void) fprintf (pfile, "255\n");
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++)
(void) fprintf (pfile, "%.9e %.9e %.9e\n",
*(picturer + ip * 600 + jp),
*(pictureg + ip * 600 + jp),
*(pictureb + ip * 600 + jp));
(void) fclose (pfile);
(void) printf ("Info: %s wrote file %s\n", *(argv + 0), &buf[0]);
pixc1old = pixc0old = 0.0;
for (inormalization = 0; inormalization < 100; inormalization++) {
vmax = 0.0;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
vmax = vmax < *(picturer + ip * 600 + jp) ? *(picturer + ip * 600 + jp) : vmax;
vmax = vmax < *(pictureg + ip * 600 + jp) ? *(pictureg + ip * 600 + jp) : vmax;
vmax = vmax < *(pictureb + ip * 600 + jp) ? *(pictureb + ip * 600 + jp) : vmax;
}
vmin = vmax;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
vmin = vmin > *(picturer + ip * 600 + jp) ? *(picturer + ip * 600 + jp) : vmin;
vmin = vmin > *(pictureg + ip * 600 + jp) ? *(pictureg + ip * 600 + jp) : vmin;
vmin = vmin > *(pictureb + ip * 600 + jp) ? *(pictureb + ip * 600 + jp) : vmin;
}
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
*(picturer + ip * 600 + jp) -= vmin;
*(pictureg + ip * 600 + jp) -= vmin;
*(pictureb + ip * 600 + jp) -= vmin;
}
vmax -= vmin;
picscale = 256.0 / vmax;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
*(picturer + ip * 600 + jp) *= picscale;
*(pictureg + ip * 600 + jp) *= picscale;
*(pictureb + ip * 600 + jp) *= picscale;
}
pixstandarddeviation = pixaverage = 0.0;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
vmax = 0.0;
vmax = vmax < *(picturer + ip * 600 + jp) ? *(picturer + ip * 600 + jp) : vmax;
vmax = vmax < *(pictureg + ip * 600 + jp) ? *(pictureg + ip * 600 + jp) : vmax;
vmax = vmax < *(pictureb + ip * 600 + jp) ? *(pictureb + ip * 600 + jp) : vmax;
pixstandarddeviation += vmax * vmax;
pixaverage += vmax;
}
pixstandarddeviation -= pixaverage * pixaverage / (double) ((irowbeyond - irowmin) * (jcolbeyond - jcolmin));
pixstandarddeviation /= (double) ((irowbeyond - irowmin) * (jcolbeyond - jcolmin) - 1);
pixstandarddeviation = pixstandarddeviation < (double) 1.0e-300 ? 0.0 : pixstandarddeviation;
pixstandarddeviation = sqrt (pixstandarddeviation);
pixaverage /= (double) ((irowbeyond - irowmin) * (jcolbeyond - jcolmin));
pixdeviation = pixstandarddeviation * 3.0;
if (pixdeviation < (double) 1.0e-300) {
pixc1 = 0.0;
pixc0 = 128.0;
}
else {
pixc1 = 128.0 / pixdeviation;
pixc0 = (pixdeviation - pixaverage) * pixc1;
}
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
pixvalmapped = *(picturer + ip * 600 + jp) * pixc1 + pixc0;
if (pixvalmapped < 0.0)
pixvalmapped = 0.0;
else if (pixvalmapped > 256.0)
pixvalmapped = 256.0;
*(picturer + ip * 600 + jp) = pixvalmapped;
pixvalmapped = *(pictureg + ip * 600 + jp) * pixc1 + pixc0;
if (pixvalmapped < 0.0)
pixvalmapped = 0.0;
else if (pixvalmapped > 256.0)
pixvalmapped = 256.0;
*(pictureg + ip * 600 + jp) = pixvalmapped;
pixvalmapped = *(pictureb + ip * 600 + jp) * pixc1 + pixc0;
if (pixvalmapped < 0.0)
pixvalmapped = 0.0;
else if (pixvalmapped > 256.0)
pixvalmapped = 256.0;
*(pictureb + ip * 600 + jp) = pixvalmapped;
}
epsiter = sqrt ((pixc1 - pixc1old) * (pixc1 - pixc1old) + (pixc0 - pixc0old) * (pixc0 - pixc0old));
pixc1old = pixc1;
pixc0old = pixc0;
if (epsiter < (double) 1.0e-12)
break;
}
(void) printf ("Info: %s performed %u normalization cycles\n", *(argv + 0), inormalization);
vmax = 0.0;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
vmax = vmax < *(picturer + ip * 600 + jp) ? *(picturer + ip * 600 + jp) : vmax;
vmax = vmax < *(pictureg + ip * 600 + jp) ? *(pictureg + ip * 600 + jp) : vmax;
vmax = vmax < *(pictureb + ip * 600 + jp) ? *(pictureb + ip * 600 + jp) : vmax;
}
vmin = vmax;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
vmin = vmin > *(picturer + ip * 600 + jp) ? *(picturer + ip * 600 + jp) : vmin;
vmin = vmin > *(pictureg + ip * 600 + jp) ? *(pictureg + ip * 600 + jp) : vmin;
vmin = vmin > *(pictureb + ip * 600 + jp) ? *(pictureb + ip * 600 + jp) : vmin;
}
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
*(picturer + ip * 600 + jp) -= vmin;
*(pictureg + ip * 600 + jp) -= vmin;
*(pictureb + ip * 600 + jp) -= vmin;
}
vmax -= vmin;
picscale = 256.0 / vmax;
for (ip = irowmin; ip < irowbeyond; ip++)
for (jp = jcolmin; jp < jcolbeyond; jp++) {
*(picturer + ip * 600 + jp) *= picscale;
*(pictureg + ip * 600 + jp) *= picscale;
*(pictureb + ip * 600 + jp) *= picscale;
}
(void) sprintf (&buf[0], "DiffractionI%uD%uW%02uS%03uB%02uR%05uV%03uL%03uH%03uL%03uOverview.ppm",
idiaphragmpos, isrcdisplacement, kwavelength, nsrcpoints, nblades,
nrad, irowmin, irowbeyond, jcolmin, jcolbeyond);
pfile = fopen (&buf[0], "w");
if (pfile == (FILE *) 0) {
(void) printf ("Unable to open %s.\n", &buf[0]);
return (1);
}
(void) fprintf (pfile, "P3\n");
(void) fprintf (pfile, "600 400\n");
(void) fprintf (pfile, "255\n");
for (ip = 0; ip < 400; ip++)
for (jp = 0; jp < 600; jp++) {
apixr = *(picturer + ip * 600 + jp);
apixg = *(pictureg + ip * 600 + jp);
apixb = *(pictureb + ip * 600 + jp);
apixr = apixr > 256.0 ? 256.0 : apixr;
apixg = apixg > 256.0 ? 256.0 : apixg;
apixb = apixb > 256.0 ? 256.0 : apixb;
pixr = (unsigned int) apixr;
pixg = (unsigned int) apixg;
pixb = (unsigned int) apixb;
pixr = pixr > 255 ? 255 : pixr;
pixg = pixg > 255 ? 255 : pixg;
pixb = pixb > 255 ? 255 : pixb;
(void) fprintf (pfile, "%u %u %u\n", pixr, pixg, pixb);
}
(void) fclose (pfile);
(void) printf ("Info: %s wrote file %s\n", *(argv + 0), &buf[0]);
}
(void) free (pictureb);
(void) free (pictureg);
(void) free (picturer);
#if (defined (MYUSEMPI))
(void) MPI_Finalize ();
#endif
return (0);
}
static void
heapsortsrcpoints (double *pdist, unsigned int *pirow, unsigned int *picol, unsigned int n)
{
unsigned int i, j, irowaux, icolaux, k;
double distaux;
for (i = n >> 1; i > 0; i--) {
j = i - 1;
distaux = *(pdist + j);
irowaux = *(pirow + j);
icolaux = *(picol + j);
k = (j << 1) + 1;
while (k < n) {
if (k + 1 < n)
if (*(pdist + k + 1) > *(pdist + k + 0))
k++;
if (distaux >= *(pdist + k))
break;
*(pdist + j) = *(pdist + k);
*(pirow + j) = *(pirow + k);
*(picol + j) = *(picol + k);
j = k;
k = (j << 1) + 1;
}
*(pdist + j) = distaux;
*(pirow + j) = irowaux;
*(picol + j) = icolaux;
}
for (i = n; i > 1; i--) {
distaux = *(pdist + i - 1);
*(pdist + i - 1) = *(pdist + 0);
irowaux = *(pirow + i - 1);
*(pirow + i - 1) = *(pirow + 0);
icolaux = *(picol + i - 1);
*(picol + i - 1) = *(picol + 0);
j = 0;
k = 1;
while (k + 1 < i) {
if (k + 2 < i)
if (*(pdist + k + 1) > *(pdist + k + 0))
k++;
if (distaux >= *(pdist + k))
break;
*(pdist + j) = *(pdist + k);
*(pirow + j) = *(pirow + k);
*(picol + j) = *(picol + k);
j = k;
k = (j << 1) + 1;
}
*(pdist + j) = distaux;
*(pirow + j) = irowaux;
*(picol + j) = icolaux;
}
return;
}
static double
x2nf (const char *str)
{
unsigned int u;
double y;
if (sscanf (str, "%x", &u) != 1) {
(void) printf ("Unrecoverable error.\n");
(void) exit (1);
}
if (u > 255) {
(void) printf ("Unrecoverable error.\n");
(void) exit (1);
}
y = (double) u / 255.0;
return (y);
}
Stephan K.H. Seidl