#!/usr/bin/perl -w  

#####################################################
# Program to synthesize astronomical spectra, and convert them into
# waveforms, ready for conversion into
# sound. Uses the Fast Fourier Transform (Math::FFT) Perl module available from CPAN.
#
# run as follows:
# ./playlines.pl <lines.d  > spec.dat
# ./sox spec.dat spec.wav
#
# 
# The synthesized spectrum can consist of any number of emission lines,
# and/or of black body (Planck law) or power-law continuum emission.
#
# Emission-line parameters are read from a text file. Continuum parameters
# are set by editing this program.
#
# Parameters you might want to edit are enclosed in ###############
#
# Paul Francis, ANU, 7th Sept 2005.
#####################################################

use Math::FFT;
my $PI = 3.1415926539;

#####################################################
# Length of time (in 11.88 sec blocks) of sample clip. Integer
# If you want a longer sound clip, increase lclip.

my $lclip = 1;
#####################################################

# FFT size. Must be 2**n.
my $N = 524288;
 
#####################################################  
# Frequency of H-alpha (Hz): 261.63 = middle C.
# This sets the conversion of light frequencies into sound frequencies. It is
# set up so that the H-alpha emission line (656.3nm wavelength) sounds like
# middle C.

my $freq_alpha = 261.63;  
#####################################################
  
# Set up the power spectrum. 
  
my ($series);
for (my $k=0; $k<$N; $k++) {
    $spec->[$k] = 0.0;
    $series->[2*$k] = 0.0;
    $series->[2*$k+1] = 0.0;
}

# Read in list of emission lines. Wavelength, peak flux and velocity width.
# You need to prepare a text file containing these values, one per line, with no
# carriage return at the end of the last line. If you don't want emission lines,
# use a file with only one line and give it zero peak flux.

my $nlines = 0;
while(<STDIN>){
   my @numbers = split;
   $wave[$nlines] = $numbers[0];
   $freq[$nlines] = $freq_alpha*6563.0/($wave[$nlines]); # convert to Hz
   $flux[$nlines] = $numbers[1];
   $amp[$nlines] = 2.32E-8*$flux[$nlines]*$wave[$nlines]*$wave[$nlines]; # Convert to Fnu
   $vwidth[$nlines] = $numbers[2];
   $nlines++;
}

# Add Gaussian lines

for (my $j=0; $j < $nlines; $j++){

  $centre = int(11.88*$freq[$j]); # Convert from Hz to bin. Centre on bin
  my $width = 1.20*$centre*$vwidth[$j]/300000.0; # Convert to frequency FWHM.
  my $peak = rand(0.8)+0.1;
  for (my $k=0; $k<$N; $k++) {
      my $lo = $centre - 5*$width;
      my $hi = $centre + 5*$width;
      if ($k < $hi && $k > $lo){
        $spec->[$k] = $amp[$j]*exp(-1.0*(($k-$centre)/$width)**2); # Gaussian 
      }
  }
}


# Now add the continuum.
# You can add a power-law or black body continuum. If you don't want these, set the
# amplitudes to sero.


# Power-law continuum 

##################################
my $poweramp = 0.0;
##################################

for (my $k=1; $k<$N; $k++) {
       $spec->[$k] += $poweramp*$k**(-0.7); 
}

# Planck Law continuum.

##################################
# Set the temperature T (in Kelvin) and amplitude $plamp.

my $plamp = 0.0;
my $T=2200.0;
##################################


my $pmax = 0.0;

for (my $k=0; $k<$N; $k++) {

    $planck[$k] = ($k**3)/(exp(39.2*200.0*$k/($T*$freq_alpha)-1));
    
    if ($planck[$k] > $pmax){
       $pmax = $planck[$k];
    }
}

for (my $k=0; $k<$N; $k++) {
    $spec->[$k] += $plamp*$planck[$k]/$pmax;
}


# Convert to amplitude spectrum and add in random phases

for (my $k=0; $k<$N; $k++) {
    my $phase = rand(2*$PI);
#    if ($k < 10000){
#       printf "%g %g \n", $k, $spec->[$k];
#    }
    $spec->[$k] = sqrt($spec->[$k]);
    $series->[2*$k] += cos($phase)*$spec->[$k];
    $series->[2*$k+1] += sin($phase)*$spec->[$k]; 
}

# Do the fourier transform

my $fft = new Math::FFT($series);
my $coeff = $fft->cdft();

# Normalise

my $max = 0;
my $min = 0;
my $norm = 1;
for (my $k=0; $k<$N; $k++) {
   if ($coeff->[2*$k] > $max){
      $max = $coeff->[2*$k];
   }
   if ($coeff->[2*$k] < $min){
      $min = $coeff->[2*$k];
   }
}
if ($max + $min > 0.0){
   $norm = $max;
} else {
   $norm = -1.0*$min;
}


# Output waveform.

print "; Sample Rate 44100\n"; # Correct value for garageband import.
my $nstart=0;
my $count = 0;
for (my $j=0; $j<$lclip; $j++){
   for (my $k=0; $k<$N; $k++) {
       my $temp = $coeff->[2*$k]/$norm;
       printf "%g %g \n", $count, $temp;    
       $count++;
   }
}