/* bb-spec-impulse.c Create an impulse with the same spectrum as the Sun's black body radiation. * * Tab = 4 * * Tweak initial states of some variables and compile with GCC: * * cc bb-spec-impulse.c -lm -o bb-spec-impulse * * "-lm" is apparently needed to link libm.a when including math.h. * * No command line parameters - just run it: * * ./bb-spec-impulse * * ./bb-spec-impulse > bb-spec-impulse.txt * * Robin Whittle Melbourne Australia 12 September 2006 rw@firstpr.com.au http://astroneu.com * * Sample time = 10^-17 seconds = 10 attoseconds = 0.01 femtoseconds = 3 nanometres. * Sample rate = 10^17 = 100 petahertz = 0.1 exahertz. * Maximum frequency (two sinewaves in output array) = 25 petahertz = 12nm wavelength. * * 16384 samples. * Total time = 1.6384 x 10^-13 seconds = 0.16384 picoseconds = 49.152 microns. * * Composed of 4095 sine waves, of which the lowest puts two full cycles in the * array. All sine waves have their zero phase point at the centre of the array. * This way we get a compact, phase coherent, impulse response. * */ #include #include #include /*--------------------------------------------------------------------------------*/ int main(int argc, char *argv[]) { /* Photosphere temperature in Kelvin. */ double tempp = 5700; /* Frequency range - will be divided into 4096 steps. * 25 petahertz = 2.5 x 10^16 Hz. */ double frange = 2.5e+16; /* Planck's constant: 6.626 x 10^-34 Joule seconds */ double cplanck = 6.626e-34; /* Speed of light: 2.99 x 10^8 metres/sec. */ double cc = 2.99e+8; /* PI and Boltzmann's Constant: 1.38 x 10^-23 Joules / Kelvin. */ double cpi = 3.141592653589793; double ck = 1.38e-23; /* Frequency (Hz) of the current Planck curve calculation. */ double freq; /* Energy at that frequency. */ double energy; /* Raw sine value, before multiplication by level. */ double sinraw; /* Value of sinewave to be written into an element * of tdout[]. */ double sinval; /* Output array - 16384 time domain samples. */ double tdout[16384]; /* Indexes into the left and right halves of * tdout[] - used when writing the sine wave. */ int pl, pr; /* What frequency step are we on now? 1 to 4095. */ int step; /* Temporary integer. */ int t1; /* Handles for output files. */ FILE *outbin; /* 16 bit integer for writing 16 bit output file. */ short int wordout; /* Open a file for binary writing. The ".pcm" * is intended to be helpful when loading the * file into an audio editing program. */ if (( outbin = fopen("outfile-bin.pcm", "wb"))==NULL) { printf("Cannot open output file!\n"); exit(1); } /* Clear the output array. */ for (t1 = 0; t1 < 16384; t1++) { tdout[t1] = 0; } for (step = 1; step < 4096; step++) { /* Set the frequency for this calculation. */ freq = ((float) step / 4096) * frange; /* Calculate energy: * * 2 * cplanck * (freq^3) * ----------------------------------------------- * c^2 * e^( (cplank * freq / (ck * temp)) - 1) * */ energy = 2.0 * cplanck * (freq*freq*freq) / ( cc * cc * pow(2.7182818284, ( (cplanck * freq / (ck * tempp)) - 1) ) ); printf("Step = %d Freq = %e Energy = %e \n", step, freq, energy); /* Should result in lines including: * * Step = 1 Freq = 6.103516e+12 Energy = 8.702521e-12 * * Step = 56 Freq = 3.417969e+14 Energy = 9.039307e-08 * Step = 57 Freq = 3.479004e+14 Energy = 9.054553e-08 * Step = 58 Freq = 3.540039e+14 Energy = 9.061452e-08 <<< Peak * Step = 59 Freq = 3.601074e+14 Energy = 9.060272e-08 * Step = 60 Freq = 3.662109e+14 Energy = 9.051287e-08 * * Step = 360 Freq = 2.197266e+15 Energy = 3.913618e-12 * * Note that the peak energy per frequency increment is * around 3.54 x 10^14 Hz. This is 844 nanometres. * * The peak wavelength (according to other sources) is around * 501 nanometres. * * The discrepency is caused by the wavelength approach * calculating energy within a particular (in principle * very small) increment of wavelength. At shorter * wavelengths this captures a greater range of frequencies * so leading to a shorter wavelength, higher frequency, peak * than if we were calculating the energy contained within * a short range of frequency. * * This is a little tricky. To make it more concrete, * let us consider a 1% increment in frequency. * * Freq: 2.99 x 10^14 3.0199 x 10^14 * Wavelength: 1000.000nm 990.099nm = 9.901nm * * Freq: 5.98 x 10^14 6.0398 x 10^14 * Wavelength: 500nm 495.049nm = 4.960nm * * */ /* Loop to write the sine wave, at the required * level = square root of energy. * * Sine wave starts at phase 0 at the centre-point * of the array and goes forwards, and backwards, * from there - with the backwards part being * of opposite polarity. * * Scale the level by some arbitrary fudge factor * so we are not getting such extremely low levels. * * Our frequency steps are organised so that at the * highest frequency: * * (4095 / 4096) * frange * * we approach the Nyquist frequency of the array - one * sample is for 0, the next for 180 degrees (PI) etc. * * The first frequency step is for frange / 4096. * So the frequency steps are in those multiples. * The sine calculation for the first frequency step * will write almost a full sine wave (2PI angle) over * the 8192 entries to the right, and likewise to the left * in tdout[], but negated and going from the centre to * the left. So in this frequency step, for each * sample, we increment the angle by (PI * 2 / 8192) * The second frequency step will write two sinewaves * per side, so we step (2 * PI * 2 / 8192) etc. * * So the angle for each sine calculation is: * * (t1 * step * PI / 4096) * * tdout[8192] remains unchanged as zero. * * We do not calculate the first sine value - it is 0. * So start with t1 set to 1. */ pr = 8193; pl = 8191; for (t1 = 1; t1 < 8192; t1++) { sinraw = sin(t1 * step * cpi / 4096 ); sinval = sqrt(energy) * 1000000 * sinraw; if ( (step == 10) ||(step == 60) ) { printf("t1 = %d pl = %d pr = %d sinraw = %f sinval = %f\n", t1, pl, pr, sinraw, sinval); } tdout[pl--] -= sinval; tdout[pr++] += sinval; } } /* Now write a 16 bit file for reading into an * audio program. This assumes littlendian * short integer = 16 bits. Audio programs can * usually load raw 16 bit mono files. */ for (t1 = 0; t1 < 16384; t1++) { printf("t1 = %d value = %f \n", t1, tdout[t1]); wordout = (int) tdout[t1]; fwrite(&wordout, 2, 1, outbin); } exit(0); }