/* Spectroscopic Toolkit version 1.66 by Pieter Suurmond, august 18, 2003. ST translates spectroscopic (i.e. quantumphysical) data to AIFF-file. Uses the binary database-file gfall.pbin, produced by the pconv-program that reads Robert Kurucz's latest data file (see subdirectory pconv). Play all frequencies of one element together, almost. Atomic.num: Suggested strength: Au = 7900, 0.92 Ag = 4700, 1.50 Fe = 2600, 0.18 Hg = 8000, 1.20 Sn = 5000, 1.12 Cu = 2900, 0.88 Pb = 8200, 1.40 Au+ = 7901, missing Ag+ = 4701, 1.20 Fe+ = 2601, 0.12 Hg+ = 8001, 1.10 Sn+ = 5001, 1.15 Cu+ = 2901, 0.20 Pb+ = 8201, 1.20 */ #include /* Standard headers. */ #include #include #include "ST_wavelets.h" /* Prototype addWAVELET(). */ #include "gfall.pbin.h" /* Defines structure of the binary database-file, */ /* same header as included in the 'pconv'-program. */ #include "ST_more.pbin.h" /* To extract some more out of the database. */ /* But also path to the database. */ #include "c00.h" /* To check prototypes. */ /*------------------------------------------------------------------------------------------------*/ short c01_massiveCloud(unsigned short selected, /* 100 * Atomic number + ionisation level. */ double strength, /* Linear overall amplitude (0.18=OK for Fe). */ double starttime, /* Starttime in seconds. */ double duration, /* 'effective duration' in seconds. */ ENGINE E, FILE* msg) { FILE *pbinDatabaseFP; pieterPacked pbinBuff; derivedData more; /* See ST.gfall.pbin.c. */ unsigned long linesRead = 0, linesUsed = 0, linesSkipped = 0; short e = 0; double audioHertz, a, dur, p, inflectStart, d; double loFreq = 9.99E99, /* AUDIO-freq. */ hiFreq = -9.99E99, loAmp = 9.99E99, /* Hi-Hertz(!), not audio-amplitude yet. */ hiAmp = -9.99E99; double durA, durB; /* Slope & offset. Lowest freq->1; highest freq->0.5. */ /*------------------------------------------------------------ Frequency transposition: */ double transpositionRatio = pow(2.0, -40.0); /* 40 octaves down. */ /*------------------------------------------------------------ Smoother frequency windowing: */ double lowestFreq = 20.0; /* (trapeziod instead of brickwall) */ double aboveLowestFreq = 1.25 * lowestFreq; /* Rise to 1 within 1 major third. */ double aboveHertz = aboveLowestFreq - lowestFreq; double nyquist = 0.5 * (double)getSamplerateENGINE(E); double belowNyquist = nyquist / 1.2; /* Fall to 0 within 1 minor third. */ double belowHertz = nyquist - belowNyquist; /*------------------------------------------------------- Open binary version of Kurucz's file. */ pbinDatabaseFP = fopen(kDbase, "rb"); if (pbinDatabaseFP) { /*------------------------------ CYCLE #1 (analysis) ------------------*/ while (fread(&pbinBuff, sizeof(pbinBuff), 1, pbinDatabaseFP)) { if (selected == pbinBuff.ELEM) { moreLineData(&pbinBuff, &more); /* Calc some more data. */ audioHertz = more.frequency * transpositionRatio; /* Exclude out-of-range-freqs from analysis. */ /* There may slip through some freqs that will be too weak in the second step! */ if ((audioHertz > lowestFreq) && (audioHertz < nyquist)) { if (audioHertz < loFreq) loFreq = audioHertz; if (audioHertz > hiFreq) hiFreq = audioHertz; if (more.Avalue < loAmp) loAmp = more.Avalue; if (more.Avalue > hiAmp) hiAmp = more.Avalue; } } } /*--------------------- Print statistics for the selected metals: ------------*/ if (msg) fprintf(msg, "-------- c01_massiveCloud() element %d: --------\n\ loAmp = %.8f; hiAmp = %.8f; RATIO=%.2f.\nloFreq = %.6f; hiFreq = %.6f; RATIO=%.2f.\n", selected, loAmp, hiAmp, hiAmp / loAmp, loFreq, hiFreq, hiFreq / loFreq); /*---------------- Linear duration compensation function: ------------------------------*/ /* d(loF) = 1 1 = a loF + b a = 0.5 / (loF - hiF) */ /* d(hiF) = 0.5 0.5 = a hiF + b - 0.5 = (0.5 / (loF - hiF)) hiF + b */ /* d(F) = a F + b ------------------- b = 0.5 - (0.5 hiF / (loF - hiF)) */ /* 0.5 = a (loF - hiF) */ d = loFreq - hiFreq; durA = 0.5 / d; /* A=slope and B=offset. */ durB = 0.5 - (0.5 * hiFreq / d); /*---------------- Linear bass-boost-filter-function (per element): --------------------*/ /* SAME "d" as with duration may be used: highest freq (per element) gets -6.021 dB. */ /*---------------------------- CYCLE #2: Post (audio-)wavelets into linked list. -------*/ rewind(pbinDatabaseFP); if (msg) { fprintf(msg, "1/transpositionRatio = %.1f.\n", 1.0 / transpositionRatio); fprintf(msg, "'effective duration' = %.6f s.\n", duration); } while (fread(&pbinBuff, sizeof(pbinBuff), 1, pbinDatabaseFP)) { /*----------------------------------------------- ELEMENT / ION - selector: */ if (selected == pbinBuff.ELEM) { moreLineData(&pbinBuff, &more); /* Calc some more data. */ /*----------------------------------------------------------------------*/ audioHertz = more.frequency * transpositionRatio; if ((audioHertz > lowestFreq) && (audioHertz < nyquist)) { /*------------------ Less brickwall-shaped filtering: --------------*/ if (audioHertz < aboveLowestFreq) a = (audioHertz - lowestFreq) / aboveHertz; /* Rise to 1 within 1 minor third. */ else if (audioHertz > belowNyquist) a = (nyquist - audioHertz) / belowHertz; /* Fall to 0 within 1 major third. */ else a = 1.0; /*---------------- Frequency-dependant Compensation per element: ---------------*/ /* Only lowest freq of selected element gets unit-duration and unit-amplit., */ /* all higher frequencies sound shorter. loFreq -> 1, hiFreq -> 0.5. */ d = (durA * audioHertz) + durB; /*------------------ Bass boost-filter (per element): --------------*/ a *= d; /* Highest frequency per element -6.021 dB. */ /* ----------------- Amplitude: ------------------------------------*/ /* Compress amplitude-range by approx 4th-power-root. */ a *= strength * pow(1.90e-13 * more.Avalue, 0.26); /*-------------------- Amplitude treshold (23 bit) -----------------*/ if (a >= 1.192093E-7) { /*------------------------ Panning: ---------------------------------------*/ if (pbinBuff.LOW.J == pbinBuff.UPP.J) p = 0.0; else if (pbinBuff.LOW.J > pbinBuff.UPP.J) /* Depend upon quantum-states. */ p = +0.5; else p = -0.5; if (linesUsed & 1L) /* Alternate between left and right channel. */ p += 0.118; else p -= 0.118; /*---------------- Duration & starttime: ---------------------------*/ dur = duration * d; /* Only lowest freq gets unit-duration. */ inflectStart = starttime; /* 60+: 0 to 180. */ if (addWAVELET(E, msg, /* Hz */ audioHertz, /* amp */ a, /* pan */ p, /* top */ inflectStart + (0.5 * dur), /* All start-inflections together. */ /* dur */ dur, /* Top at inflectstart + DUR/2. */ /* env */ kENV_Gaussian, /* syn */ 1)) /* 1 = sync with other oscillators. */ { if (msg) fprintf(msg, "Illegal wavelet!\n"); e = 3; goto clFile; } linesUsed++; } else { linesSkipped++; if (msg) fprintf(msg, "SKIPPED-AMP %d: Hz=%12.6f lggf=%5.2f.\n", pbinBuff.ELEM, audioHertz, pbinBuff.LOGGF); } } else { linesSkipped++; if (msg) fprintf(msg, "SKIPPED-FRQ %d: Hz=%12.6f lggf=%5.2f.\n", pbinBuff.ELEM, audioHertz, pbinBuff.LOGGF); } printLineEssential(&pbinBuff, msg); /* NULL, stdout or some FILE*. */ printLineDerived(&more, msg); } linesRead++; } if (msg) fprintf(msg, " linesUsed = %ld.\nlinesSkipped = %ld.\n", linesUsed, linesSkipped); if (linesRead != kKuruczLinesTotal) { if (msg) fprintf(msg, "Error reading binary database lines!\n"); e = 2; } clFile: fclose(pbinDatabaseFP); } else { if (msg) fprintf(msg, "Could not open binary database-file \"%s\"!\n", kDbase); e = 1; } return(e); } /*------------------------------------------------------------------------------------------------*/ short c01_shapedCloud(unsigned short selected, /* 100 * Atomic number + ionisation level. */ double strength, /* Linear overall amplitude (0.18=OK for Fe). */ double starttime, /* Starttime in seconds. */ double duration, /* 'effective duration' in seconds. */ const double* shapeFreqs, /* Array with frequencies (- earlier, + later) */ ENGINE E, FILE* msg) { FILE *pbinDatabaseFP; pieterPacked pbinBuff; derivedData more; /* See ST.gfall.pbin.c. */ unsigned long linesRead = 0, linesUsed = 0, linesSkipped = 0; short e = 0; double audioHertz, a, dur, p, inflectStart, d; double loFreq = 9.99E99, /* AUDIO-freq. */ hiFreq = -9.99E99, loAmp = 9.99E99, /* Hi-Hertz(!), not audio-amplitude yet. */ hiAmp = -9.99E99; double durA, durB; /* Slope & offset. Lowest freq->1; highest freq->0.5. */ double shF, shD, corr, shCorr; const double *shapeF, shapeRatio = 1.12; /* About a whole tone around center-frequencies. */ /*------------------------------------------------------------ Frequency transposition: */ double transpositionRatio = pow(2.0, -40.0); /* 40 octaves down. */ /*------------------------------------------------------------ Smoother frequency windowing: */ double lowestFreq = 20.0; /* (trapeziod instead of brickwall) */ double aboveLowestFreq = 1.25 * lowestFreq; /* Rise to 1 within 1 major third. */ double aboveHertz = aboveLowestFreq - lowestFreq; double nyquist = 0.5 * (double)getSamplerateENGINE(E); double belowNyquist = nyquist / 1.2; /* Fall to 0 within 1 minor third. */ double belowHertz = nyquist - belowNyquist; /*------------------------------------------------------- Open binary version of Kurucz's file. */ pbinDatabaseFP = fopen(kDbase, "rb"); if (pbinDatabaseFP) { /*------------------------------ CYCLE #1 (analysis) ------------------*/ while (fread(&pbinBuff, sizeof(pbinBuff), 1, pbinDatabaseFP)) { if (selected == pbinBuff.ELEM) { moreLineData(&pbinBuff, &more); /* Calc some more data. */ audioHertz = more.frequency * transpositionRatio; /* Exclude out-of-range-freqs from analysis. */ /* There may slip through some freqs that will be too weak in the second step! */ if ((audioHertz > lowestFreq) && (audioHertz < nyquist)) { if (audioHertz < loFreq) loFreq = audioHertz; if (audioHertz > hiFreq) hiFreq = audioHertz; if (more.Avalue < loAmp) loAmp = more.Avalue; if (more.Avalue > hiAmp) hiAmp = more.Avalue; } } } /*--------------------- Print statistics for the selected metals: ------------*/ if (msg) fprintf(msg, "-------- c01_massiveCloud() element %d: --------\n\ loAmp = %.8f; hiAmp = %.8f; RATIO=%.2f.\nloFreq = %.6f; hiFreq = %.6f; RATIO=%.2f.\n", selected, loAmp, hiAmp, hiAmp / loAmp, loFreq, hiFreq, hiFreq / loFreq); /*---------------- Linear duration compensation function: ------------------------------*/ /* d(loF) = 1 1 = a loF + b a = 0.5 / (loF - hiF) */ /* d(hiF) = 0.5 0.5 = a hiF + b - 0.5 = (0.5 / (loF - hiF)) hiF + b */ /* d(F) = a F + b ------------------- b = 0.5 - (0.5 hiF / (loF - hiF)) */ /* 0.5 = a (loF - hiF) */ d = loFreq - hiFreq; durA = 0.5 / d; /* A=slope and B=offset. */ durB = 0.5 - (0.5 * hiFreq / d); /*---------------- Linear bass-boost-filter-function (per element): --------------------*/ /* SAME "d" as with duration may be used: highest freq (per element) gets -6.021 dB. */ /*---------------------------- CYCLE #2: Post (audio-)wavelets into linked list. -------*/ rewind(pbinDatabaseFP); if (msg) { fprintf(msg, "1/transpositionRatio = %.1f.\n", 1.0 / transpositionRatio); fprintf(msg, "'effective duration' = %.6f s.\n", duration); } while (fread(&pbinBuff, sizeof(pbinBuff), 1, pbinDatabaseFP)) { /*----------------------------------------------- ELEMENT / ION - selector: */ if (selected == pbinBuff.ELEM) { moreLineData(&pbinBuff, &more); /* Calc some more data. */ /*----------------------------------------------------------------------*/ audioHertz = more.frequency * transpositionRatio; if ((audioHertz > lowestFreq) && (audioHertz < nyquist)) { /*------------------ Less brickwall-shaped filtering: --------------*/ if (audioHertz < aboveLowestFreq) a = (audioHertz - lowestFreq) / aboveHertz; /* Rise to 1 within 1 minor third. */ else if (audioHertz > belowNyquist) a = (nyquist - audioHertz) / belowHertz; /* Fall to 0 within 1 major third. */ else a = 1.0; /*---------------- Frequency-dependant Compensation per element: ---------------*/ /* Only lowest freq of selected element gets unit-duration and unit-amplit., */ /* all higher frequencies sound shorter. loFreq -> 1, hiFreq -> 0.5. */ d = (durA * audioHertz) + durB; /*------------------ Bass boost-filter (per element): --------------*/ a *= d; /* Highest frequency per element -6.021 dB. */ /* ----------------- Amplitude: ------------------------------------*/ /* Compress amplitude-range by approx 4th-power-root. */ a *= strength * pow(1.90e-13 * more.Avalue, 0.26); /*-------------------- Amplitude treshold (23 bit) -----------------*/ if (a >= 1.192093E-7) { /*------------------------ Panning: ---------------------------------------*/ if (pbinBuff.LOW.J == pbinBuff.UPP.J) p = 0.0; else if (pbinBuff.LOW.J > pbinBuff.UPP.J) /* Depend upon quantum-states. */ p = +0.5; else p = -0.5; if (linesUsed & 1L) /* Alternate between left and right channel. */ p += 0.118; else p -= 0.118; /*---------------- Shaping: -------------------------------------*/ shCorr = 0.0; /* Regions may overlap. */ shapeF = shapeFreqs; while ((shF = fabs(*shapeF)) > 1.0) /* <= 1 terminates array. */ { shD = audioHertz / shF; if (shD < 1.0) shD = 1.0 / shD; if (shD < shapeRatio) /* Within a semitone. */ { /* shD: 1.0 ... 1.059. */ /* corr = 1.0 - (log(shD) / log(shapeRatio)); 1.0 ... 0.000. */ corr = 1.0 - ((shD - 1.0) / (shapeRatio - 1.0)); corr *= corr; /* 4th power makes better shape. */ corr *= corr; if (*shapeF < 0.0) shCorr -= corr; /* Total correction may exceed 1.0 ! */ else shCorr += corr; } shapeF++; } if (msg && (shCorr != 0.0)) fprintf(msg, "shCorr = %.4f\n", shCorr); /*---------------- Duration & starttime: ---------------------------*/ dur = duration * d; /* Only lowest freq gets unit-duration. */ inflectStart = starttime + (shCorr * dur * 0.36); /* 0.60 is a bit too much. */ if (addWAVELET(E, msg, /* Hz */ audioHertz, /* amp */ a, /* pan */ p, /* top */ inflectStart + (0.5 * dur), /* All start-inflections together. */ /* dur */ dur, /* Top at inflectstart + DUR/2. */ /* env */ kENV_Gaussian, /* syn */ 1)) /* 1 = sync with other oscillators. */ { if (msg) fprintf(msg, "Illegal wavelet!\n"); e = 3; goto clFile; } linesUsed++; } else { linesSkipped++; if (msg) fprintf(msg, "SKIPPED-AMP %d: Hz=%12.6f lggf=%5.2f.\n", pbinBuff.ELEM, audioHertz, pbinBuff.LOGGF); } } else { linesSkipped++; if (msg) fprintf(msg, "SKIPPED-FRQ %d: Hz=%12.6f lggf=%5.2f.\n", pbinBuff.ELEM, audioHertz, pbinBuff.LOGGF); } printLineEssential(&pbinBuff, msg); /* NULL, stdout or some FILE*. */ printLineDerived(&more, msg); } linesRead++; } if (msg) fprintf(msg, " linesUsed = %ld.\nlinesSkipped = %ld.\n", linesUsed, linesSkipped); if (linesRead != kKuruczLinesTotal) { if (msg) fprintf(msg, "Error reading binary database lines!\n"); e = 2; } clFile: fclose(pbinDatabaseFP); } else { if (msg) fprintf(msg, "Could not open binary database-file \"%s\"!\n", kDbase); e = 1; } return(e); }