#!/usr/bin/perl

# This script computes the gain (i.e, arcseconds in CSMM space per LSB, or
# fine DAC count) of the MIPS Cryogenic Scan Mirror Mechanism (CSMM) as a
# function of CSMM position, as measured by the coarse DAC setting.  The
# data used come from a spreadsheet generated by Jerry Heim at Ball
# Aerospace.  The script uses as input a list of scan mirror deflections
# (measured with a theodolite) vs. applied voltage, which is converted to
# coarse DAC counts using the hardstop settings.  The derivative of this
# list of deflections is the CSMM gain.  A parabola is fit to the derivative
# to derive a function that describes the CSMM gain vs. coarse DAC position
# (as set by CESCANABS, for example).

# C. W. Engelbracht (2003-03-09)


use PGPLOT;
use PDL;
use PDL::Fit::Polynomial;
use PDL::Graphics::PGPLOT;


# set parameters

$file = "flight_CSMM_Test_Data.csv";
$hardstop1 = 610;
$hardstop2 = 3486;


# read data

open (FILE, $file) or die "Cannot open $file: $!";

for ($i = 0; $i < 5; $i++) {<FILE>;}

$i = 0;

while (<FILE>) {
   @line = split (',', $_);
   ($volt1[$i], $angle1[$i], $volt2[$i], $angle2[$i]) = 
      ($line[0], -$line[2], $line[6], -$line[8]);
   $i++;
}

close (FILE);

until ($volt1[$#volt1] ne "") {pop @volt1; pop @angle1;}
until ($volt2[0] ne "") {shift @volt2; shift @angle2;}


# ignore first point on secondary side because it doesn't fit well

shift @volt2; shift @angle2;


# convert volts to DAC counts

$convert = ($hardstop2 - $hardstop1 ) / ($volt1[$#volt1] - $volt1[0]);
for ($i = 0; $i <= $#volt1; $i++) {
   $count1[$i] = ($volt1[$i] - $volt1[0]) * $convert + $hardstop1;
}

$convert = ($hardstop2 - $hardstop1 ) / ($volt2[$#volt2] - $volt2[0]);
for ($i = 0; $i <= $#volt2; $i++) {
   $count2[$i] = ($volt2[$i] - $volt2[0]) * $convert + $hardstop1;
}


# extract linear part of angle vs. DAC data

$delta = 1100;
$dac1 = $hardstop1 + $delta;
$dac2 = $hardstop2 - $delta;

$j = 0;
for ($i = 0; $i <= $#count1; $i++) {
   if (($count1[$i] > $dac1) && ($count1[$i] < $dac2)) {
      $smallcount1[$j] = $count1[$i];
      $smallangle1[$j] = $angle1[$i];
      $j++;
   }
}

$j = 0;
for ($i = 0; $i <= $#count2; $i++) {
   if (($count2[$i] > $dac1) && ($count2[$i] < $dac2)) {
      $smallcount2[$j] = $count2[$i];
      $smallangle2[$j] = $angle2[$i];
      $j++;
   }
}


# fit lines to linear part of angle vs. DAC data

@fit1 = fitline (\@smallcount1, \@smallangle1);
@fit2 = fitline (\@smallcount2, \@smallangle2);
printf "Average CSMM gain fit over CESCANABS %0d to %0d:\n", $dac1, $dac2;
printf "   CE1:  %0.2f +/- %0.3f\n", $fit1[0] * 3600, $fit1[3] * 3600;
printf "   CE2:  %0.2f +/- %0.3f\n", $fit2[0] * 3600, $fit2[3] * 3600;


# compute derivatives

$deriv1[0] = ($angle1[1] - $angle1[0]) / ($count1[1] - $count1[0]) * 3600;
for ($i = 1; $i < $#count1; $i++) {
   $deriv1[$i] = ($angle1[$i+1] - $angle1[$i-1])
      / ($count1[$i+1] - $count1[$i-1]) * 3600;
}
$deriv1[$#count1] = ($angle1[$#count1] - $angle1[$#count1-1])
   / ($count1[$#count1] - $count1[$#count1-1]) * 3600;

$deriv2[0] = ($angle2[1] - $angle2[0]) / ($count2[1] - $count2[0]) * 3600;
for ($i = 1; $i < $#count2; $i++) {
   $deriv2[$i] = ($angle2[$i+1] - $angle2[$i-1])
      / ($count2[$i+1] - $count2[$i-1]) * 3600;
}
$deriv2[$#count2] = ($angle2[$#count2] - $angle2[$#count2-1])
   / ($count2[$#count2] - $count2[$#count2-1]) * 3600;


# fit parabola to derivatives

$count1 = pdl (@count1);
$deriv1 = pdl (@deriv1);
($gainfit1, $gaincoeffs1) = fitpoly1d ($count1, $deriv1, 3);
@gainfit1 = list $gainfit1;

printf "CE1 gain is %0.3e + %0.3ex + %0.3ex^2\n", at ($gaincoeffs1, (0)),
   at ($gaincoeffs1, (1)), at ($gaincoeffs1, (2));

$count2 = pdl (@count2);
$deriv2 = pdl (@deriv2);
($gainfit2, $gaincoeffs2) = fitpoly1d ($count2, $deriv2, 3);
@gainfit2 = list $gainfit2;

printf "CE2 gain is %0.3e + %0.3ex + %0.3ex^2\n", at ($gaincoeffs2, (0)),
   at ($gaincoeffs2, (1)), at ($gaincoeffs2, (2));
print "where x is the SCANABS value\n";

printf "e.g., CE1 gain at SCANABS=2007 is %0.2f\n", at ($gaincoeffs1, (0)) +
   at ($gaincoeffs1, (1)) * 2007 + at ($gaincoeffs1, (2)) * 2007 ** 2;


# renormalize angle measurements so they're centered on 0

$dac_center = ($hardstop2 + $hardstop1) / 2;
$angle1_center = $fit1[0] * $dac_center + $fit1[1];
$angle2_center = $fit2[0] * $dac_center + $fit2[1];
for ($i = 0; $i <= $#angle1; $i++) {$angle1[$i] -= $angle1_center;}
for ($i = 0; $i <= $#angle2; $i++) {$angle2[$i] -= $angle2_center;}
$fit1[1] -= $angle1_center;
$fit2[1] -= $angle2_center;


# plot data

pgopen ('?');
pgqinf ('type', $device, $junk);
if ($device ne "XSERVE") {pgpap (6, 1);}
pgsubp (2, 2);
pgsch (2);

($cmin, $cmax) = minmax (@count1);
($amin, $amax) = minmax (@angle1);
pgenv ($cmin, $cmax, $amin, $amax, 0, 0);
pglab ("coarse DAC (counts)", "CSMM angle (degrees)", "CE 1");
pgsci (2);
pgpt ($#count1+1, \@count1, \@angle1, 1);
pgsci (4);
pgline (2, [$cmin,$cmax], [$fit1[0]*$cmin+$fit1[1],$fit1[0]*$cmax+$fit1[1]]);
pgsci (1);

($cmin, $cmax) = minmax (@count2);
($amin, $amax) = minmax (@angle2);
pgenv ($cmin, $cmax, $amin, $amax, 0, 0);
pglab ("coarse DAC (counts)", "CSMM angle (degrees)", "CE 2");
pgsci (2);
pgpt ($#count2+1, \@count2, \@angle2, 1);
pgsci (4);
pgline (2, [$cmin,$cmax], [$fit2[0]*$cmin+$fit2[1],$fit2[0]*$cmax+$fit2[1]]);
pgsci (1);

($cmin, $cmax) = minmax (@count1);
($dmin, $dmax) = minmax (@deriv1);
pgenv ($cmin, $cmax, $dmin, $dmax, 0, 0);
pglab ("coarse DAC (counts)", "CSMM Gain (\"/LSB)", "");
pgsci (2);
pgpt ($#count1+1, \@count1, \@deriv1, 1);
pgsci (4);
pgline ($#count1+1, \@count1, \@gainfit1);
pgsci (1);

($cmin, $cmax) = minmax (@count2);
($dmin, $dmax) = minmax (@deriv2);
pgenv ($cmin, $cmax, $dmin, $dmax, 0, 0);
pglab ("coarse DAC (counts)", "CSMM Gain (\"/LSB)", "");
pgsci (2);
pgpt ($#count2+1, \@count2, \@deriv2, 1);
pgsci (4);
pgline ($#count2+1, \@count2, \@gainfit2);
pgsci (1);

pgend;




# SUBROUTINES

# compute minimum and maximum values of an array

sub minmax {
   my ($min, $max, $i);
   my @array = @_;
   $min = $array[0];
   $max = $array[0];
   for ($i = 1; $i <= $#array; $i++) {
      if ($array[$i] < $min) {$min = $array[$i];}
      if ($array[$i] > $max) {$max = $array[$i];}
   }
   return ($min, $max);
}


# fit a line to input x and y arrays - return slope, intercept, and errors

sub fitline {
   my @x = @{$_[0]}; my @y = @{$_[1]};
   my $sumx2; my $sumx; my $sumxy; my $sumy; my $sumdev2;
   my $i; my $delta; my $slope; my $intercept;
   my $sigma_slope; my $sigma_points;
   for ($i = 0; $i <= $#x; $i++) {
      $sumx += $x[$i];
      $sumx2 += $x[$i]**2;
      $sumxy += $x[$i] * $y[$i];
      $sumy += $y[$i];
   }
   $delta = ($#y+1) * $sumx2 - $sumx**2;
   if ($delta == 0) {return (0, 0, 0, 0);}
   $slope = (($#y+1) * $sumxy - $sumx * $sumy) / $delta;
   $intercept = ($sumx2 * $sumy - $sumx * $sumxy) / $delta;
   for ($i = 0; $i <= $#x; $i++) {
      $sumdev2 += ($y[$i] - $intercept - $slope * $x[$i]) ** 2;
   }
   $sigma_points = sqrt (1 / ($#y-1) * $sumdev2);
   $sigma_slope = sqrt (($#y+1) * $sigma_points**2 / $delta);
   return ($slope, $intercept, $sigma_points, $sigma_slope);
}