#!/usr/bin/perl
# 
# generate some fake-but-reasonable-looking data which shows
#   the received power of radio signals from a set of 5 stations
#   as a function of Q values.
#
# MWR 1/30/2020

$debug = 1;

# avoid dividing by zero
$TINY = 1.0E-10;
$BIG = 1.0/$TINY;


$mypi = 3.14159;

# range of frequencies over which to compute power (MHz)
$min_freq = 40.0;
$max_freq = 140.0;
$d_freq = 1.0;
# station frequencies (MHz) -- we are interested in the central one
@channels = (90.0, 90.5, 91.0, 91.5, 92.0);
$omega0 = $channels[2];

# arbitrary max amplitude
$A = 5800.0;

# first, we print out fixed values for the stations, 
#    giving each an equal weight of 1.0
$outfile = "./make_radio_data.out1";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
foreach $channel (@channels) {
  printf OUTFILE " omega %12.5f  power %12.5f \n", $channel, 1.0;
}
close(OUTFILE);




# range of frequency to use in calculations
$n_per_omega = 5500;
$n_omega = 2.0;
$start_omega = 0.5;
$end_omega = $n_omega*$omega0;
$d_omega = ($omega0/$n_per_omega);


$outfile = "./make_radio_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");

@Q_values = (1.5, 90.0, 150.0, 280.0, 800.0, 3000.0);
foreach $Q (@Q_values) {

		  printf OUTFILE "# Q = $Q  \n";
		  printf OUTFILE "# \n";

        # compute the amplitude at the central freqency
        $top = $A;
        $bot = (($omega0**2.0*$omega0**2.0)/($Q**2.0));
        $bot = sqrt($bot);
        $central_amp = $top / $bot;

		  for ($omega = $start_omega; $omega <= $end_omega; $omega += $d_omega) {

			 # amplitude
			 $top = $A;
			 $bot = (($omega0**2.0 - $omega**2.0)**2.0) + 
                              (($omega0**2.0*$omega**2.0)/($Q**2.0));
			 $bot = sqrt($bot);

			 $amp = $top / $bot;
          # normalize to value at central frequency
#          $amp /= $central_amp;

			 printf OUTFILE " omega %12.5f  top %12.5f  bot %12.5f  amp %12.5f   norm %12.5f  \n", 
						  $omega, $top, $bot, $amp, ($amp*$amp)/($central_amp*$central_amp);
			
		  }
        printf OUTFILE " \n";
        printf OUTFILE " \n";

        # and now repeat, but just for the frequencies of the stations.
        foreach $omega (@channels) {

			 $top = $A;
			 $bot = (($omega0**2.0 - $omega**2.0)**2.0) + 
                              (($omega0**2.0*$omega**2.0)/($Q**2.0));
			 $bot = sqrt($bot);

			 $amp = $top / $bot;
          # normalize to value at central frequency
#          $amp /= $central_amp;

			 printf OUTFILE " omega %12.5f  top %12.5f  bot %12.5f  amp %12.5f  norm %12.5f  \n", 
						  $omega, $top, $bot, $amp, ($amp*$amp)/($central_amp*$central_amp);
			
		  }
        printf OUTFILE " \n";
        printf OUTFILE " \n";

}



close(OUTFILE);


exit 0;

