/* 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);
}