/* 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 #include #include #include #include #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 #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], " "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " \\\n "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " \\\n "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " "); (void) strcat (&buf[0], " "); (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); }