#!/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 <> 5:
                sys.stderr.write ('usage: refl_calc.py open_comp_file short_comp_file load_comp_file in_file out_file\n')
                sys.exit (1)

	open_comp_file = args[0]
	short_comp_file = args[1]
	load_comp_file = args[2]
	in_file = args[3]
	out_file = args[4]

	open_comp_fh = open(open_comp_file, "r")
	short_comp_fh = open(short_comp_file, "r")
	load_comp_fh = open(load_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_files

	   # 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 parse a line of OPEN compensation data
	   open_comp_line = open_comp_fh.readline()
           field = open_comp_line[0:open_comp_line.find(sep)]
           open_comp_freq = float(field)
           open_comp_line = open_comp_line[len(field)+1:len(open_comp_line)]
           field = open_comp_line[0:open_comp_line.find(sep)]
           open_comp_mag = float(field)
           open_comp_line = open_comp_line[len(field)+1:len(open_comp_line)]
           field = open_comp_line[0:open_comp_line.find(sep)]
           open_comp_phase = float(field)
           open_comp_line = open_comp_line[len(field)+1:len(open_comp_line)]
           open_comp_sphase = float(open_comp_line)

	   # read and parse a line of SHORT compensation data
           short_comp_line = short_comp_fh.readline()
           field = short_comp_line[0:short_comp_line.find(sep)]
           short_comp_freq = float(field)
           short_comp_line = short_comp_line[len(field)+1:len(short_comp_line)]
           field = short_comp_line[0:short_comp_line.find(sep)]
           short_comp_mag = float(field)
           short_comp_line = short_comp_line[len(field)+1:len(short_comp_line)]
           field = short_comp_line[0:short_comp_line.find(sep)]
           short_comp_phase = float(field)
           short_comp_line = short_comp_line[len(field)+1:len(short_comp_line)]
           short_comp_sphase = float(short_comp_line)

	   # read and parse a line of LOAD compensation data
           load_comp_line = load_comp_fh.readline()
           field = load_comp_line[0:load_comp_line.find(sep)]
           load_comp_freq = float(field)
           load_comp_line = load_comp_line[len(field)+1:len(load_comp_line)]
           field = load_comp_line[0:load_comp_line.find(sep)]
           load_comp_mag = float(field)
           load_comp_line = load_comp_line[len(field)+1:len(load_comp_line)]
           field = load_comp_line[0:load_comp_line.find(sep)]
           load_comp_phase = float(field)
           load_comp_line = load_comp_line[len(field)+1:len(load_comp_line)]
           load_comp_sphase = float(load_comp_line)


	   # check that all frequencies match
	   if (in_freq <> open_comp_freq) | (in_freq <> short_comp_freq ) | ( in_freq <> load_comp_freq ) :
	      print "FILES OUT OF SYNC! Rerun measurements"
	      in_line = ''
	   else:
	      in_point = complex((in_mag * in_phase), (in_mag * in_sphase))

	      open_comp_point = complex((open_comp_mag * open_comp_phase), (open_comp_mag * open_comp_sphase))

	      short_comp_point = complex((short_comp_mag * short_comp_phase), (short_comp_mag * short_comp_sphase))

              load_comp_point = complex((load_comp_mag * load_comp_phase), (load_comp_mag * load_comp_sphase))


	      # From McDermott: (a) = actual    (m) = measured
	      #
	      # s11(a) = ( s11(m) - Ed ) / ( Es * ( s11(m) - Ed ) + Et )
	      #
	      # Ed = s11(m)(load)
	      #
	      # Es = ( 2 * s11(m)(load) - ( s11(m)(short) + s11(m)(open) ) ) /
	      #              ( s11(m)(short) - s11(m)(open) )
	      #
	      # Et = ( 2 * ( s11(m)(open) + s11(m)(load) ) * ( s11(m)(short) + 
	      #      s11(m)(load) ) )  /  ( s11(m)(short) - s11(m)(open) )
	      #
	      # here, s11(m)(load)  is  load_comp_point
	      #       s11(m)(short) is  short_comp_point
	      #       s11(m)(open)  is  open_comp_point
	      #

	      Ed = load_comp_point
	      Es = ( 2 * load_comp_point - ( short_comp_point + open_comp_point ) ) / ( short_comp_point - open_comp_point )
	      Et = ( 2 * ( open_comp_point + load_comp_point ) * ( short_comp_point + load_comp_point ) ) / ( short_comp_point - open_comp_point )

	      dut_point = ( in_point - Ed ) / ( Es * ( in_point - Ed ) + Et )

	      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:])