#!/bin/env /usr/bin/python

# inputs: <calibration (transmission through) file> <measured DUT file>
# output: compensated DUT file
#
# in_file measured DUT data
# comp_file calibration data using straight through transmissions, the 'bullet'
# out_file  output results, measured data compensated with through data

import sys, math

def main (args):

        nargs = len (args)
        if nargs <> 3:
                sys.stderr.write ('usage: trans_calc.py comp_file in_file out_file\n')
                sys.exit (1)

	comp_file = args[0]
	in_file = args[1]
	out_file = args[2]

	comp_fh = open(comp_file, "r")
	in_fh = open(in_file, "r")
	out_fh = open(out_file, "w")

	sep = ' '
	in_line = in_fh.readline()
	while in_line <> '': # in_file can be shorter than comp_file

           # parse a line of uncompensated measured data
           field = in_line[0:in_line.find(sep)]
           in_freq = float(field)
           in_line = in_line[len(field)+1:len(in_line)]
           field = in_line[0:in_line.find(sep)]
           in_mag = float(field)
           in_line = in_line[len(field)+1:len(in_line)]
           field = in_line[0:in_line.find(sep)]
           in_phase = float(field)
           in_line = in_line[len(field)+1:len(in_line)]
           in_sphase = float(in_line)

	   # read and arse a line of compensation data
	   comp_line = comp_fh.readline()
           field = comp_line[0:comp_line.find(sep)]
           comp_freq = float(field)
           comp_line = comp_line[len(field)+1:len(comp_line)]
           field = comp_line[0:comp_line.find(sep)]
           comp_mag = float(field)
           comp_line = comp_line[len(field)+1:len(comp_line)]
           field = comp_line[0:comp_line.find(sep)]
           comp_phase = float(field)
           comp_line = comp_line[len(field)+1:len(comp_line)]
           comp_sphase = float(comp_line)


	   if in_freq <> comp_freq:
	      print "FILES OUT OF SYNC! Rerun measurements"
	      in_line = ''
	   else:
	      in_point = complex((in_mag * in_phase), (in_mag * in_sphase))
	      comp_point = complex((comp_mag * comp_phase), (comp_mag * comp_sphase))

	      # From McDermott: s21(actual) = s21(measured) / CAL21
	      dut_point = in_point / comp_point

	      dut_mag = math.sqrt(dut_point.real**2. + dut_point.imag**2.)
	      dut_phase = dut_point.real / dut_mag
	      dut_sphase = dut_point.imag / dut_mag

	      out_fh.write(str(in_freq)+' '+str(dut_mag)+' '+str(dut_phase)+' '+str(dut_sphase)+'\n')
	      in_line = in_fh.readline()


if __name__ == '__main__':
        main (sys.argv[1:])