#!/usr/bin/env python

import datetime
import cx_Oracle
from sys import argv, exit
from os import getenv, environ, popen, umask, getcwd, system
from os.path import isfile
from string import find, upper, lower, strip, split
try:
	import psycopg2
except:
	print 'psycopg2 not found.  Cannot continue.'
	exit(0)


program = 'db2cc'
author = 'Walter Brisken'
version = '0.2'
verdate = '20150730'

# Name of the Postgres VLBA database based on EVLA
databaseName = getenv("VLBAMPTS_DB")
if databaseName == None:
	print 'Error: VMBAMPTS_DB environment variable not set.  Cannot continue.'
	exit(0)

mjd0 = datetime.datetime(1858, 11, 17, 0, 0)

def usage(pgm):
	print '%s ver %s  %s  %s\n' % (program, version, author, verdate)
	print 'A program to extract cable cal values from the wideband VLBA database.\n'
	print 'Usage: %s [options] <vex file>\n' % pgm
	print 'options can be:\n'
	print '  --help'
	print '  -h         print this message and quit\n'
	print '  --verbose'
	print '  -v         print more diagnostics to screen\n'
	print '  --force'
	print '  -f         run even if means overwriting exisiting files\n'
	print 'Example:  db2cc bf115f.vex\n'
	exit(0)

# Splits a combined obs code into its proposal and segment portions
def splitobscode(exper):
	obsSeg = ''
	proposal = exper[:]

	if len(proposal) > 3:
		if proposal[0].isalpha() and proposal[1].isalpha() and proposal[2].isdigit():
			for i in range(3, len(proposal)):
				if proposal[i].isalpha():
					obsSeg = proposal[i:]
					proposal = proposal[0:i]
					break
		if proposal[0].isalpha() and proposal[1].isdigit():
			for i in range(2, len(proposal)):
				if proposal[i].isalpha():
					obsSeg = proposal[i:]
					proposal = proposal[0:i]
					break
	
	return proposal, obsSeg

# Returns obsCode, obsSeg, { ANT : [startmjd, stopmjd] }
def vexPeek(vexFile, verbose):
	
	cmd = 'vexpeek %s' % vexFile
	if verbose > 0:
		print 'Executing command: %s' % cmd

	p = popen(cmd)
	data = p.readlines()

	if len(data) == 0:
		return 'Error', 'Error', 'Error'
	
	obsCode = upper(strip(data[0]))
	obsSeg = ''

	if obsCode[0:5] == 'ERROR':
		return 'Error', 'Error', 'Error'

	if len(obsCode) > 3:
		if obsCode[0].isalpha() and obsCode[1].isalpha() and obsCode[2].isdigit():
			for i in range(3, len(obsCode)):
				if obsCode[i].isalpha():
					obsSeg = obsCode[i:]
					obsCode = obsCode[0:i]
					break;
		if obsCode[0].isalpha() and obsCode[1].isdigit():
			for i in range(2, len(obsCode)):
				if obsCode[i].isalpha():
					obsSeg = obsCode[i:]
					obsCode = obsCode[0:i]
					break;
		

	stationTimes = {}
	for d in data[1:]:
		s = split(strip(d))
		stationTimes[upper(s[0])] = [float(s[1]), float(s[2])]

	print 'This is experiment %s %s' % (obsCode, obsSeg)

	return obsCode, obsSeg, stationTimes


def mjd2vex(mjd, dateonly=False):
	md = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]
	d = int(mjd)
	s = int((mjd - d)*86400.0 + 0.5)
	dt = datetime.timedelta(d, s)
	t = mjd0 + dt
	d = t.day + md[t.month-1]
	if t.year % 4 == 0 and t.month > 2:
		d += 1
	if dateonly:
		return '%dy%03dd' % (t.year, d)
	else:
		return '%dy%03dd%02dh%02dm%02ds' % (t.year, d, t.hour, t.minute, t.second)

def mjd2date(mjd):
	mon = ['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC']
	dt = datetime.timedelta(int(mjd+0.5/86400.0), 0)
	t = mjd0 + dt
	return '%02d-%s-%04d' % (t.day, mon[t.month-1], t.year)

def Ulow(s):
	return upper(s[0]) + lower(s[1:])

def calcMJD(t):
	if type(t) == type(1.1):
		return t	# must already be an mjd
	
	# otherwise assume it is a unix time
	dt = t - mjd0
	mjd = dt.days + dt.seconds/86400.0
	return mjd

def day2mjd(year, doy):
	return (datetime.datetime(int(year), 1, 1) + datetime.timedelta(int(doy)-1) - mjd0).days

def getCableCalsFromDB(db, startMJD, stopMJD, verbose):

	start = mjd2date(startMJD)
	stop  = mjd2date(stopMJD)

	cursor = db.cursor()
	query = """
        SELECT
          split_part(hostname_timestamp, ',', 1) as hostname,
          split_part(hostname_timestamp, ',', 2)::double precision as timestamp,
          COALESCE(RTP_inphase, -999) AS RTP_inphase,
          COALESCE(RTP_quadphase, -999) AS RTP_quadphase
        FROM crosstab(
          $$ SELECT
               UPPER(SUBSTR(hostname,1,2)) || ',' || timestamp,
               monpointname,
               monpointvalue
             FROM mcdata
             WHERE devicename = 'RTPHASE' AND
                   monpointname IN ('RTP_inphase', 'RTP_quadphase') AND
                   timestamp BETWEEN %14.8f AND %14.8f
             ORDER BY 1,2 $$,
          $$ VALUES ('RTP_inphase'), ('RTP_quadphase') $$
          ) AS ct(hostname_timestamp text,
                  RTP_inphase double precision,
                  RTP_quadphase double precision)""" % (startMJD, stopMJD)

	if(verbose):
		print 'Executing query: %s' % query
	cursor.execute(query)
	ccdata = cursor.fetchall()

	return ccdata

def processVexFile(vexFile, force, verbose):

	CHALFM =  7344000.0
	CHALFP =  7346400.0
	CHALF  =  7348800.0
	CMAX   = 14692800.0

	obsCode, obsSeg, stationTimes = vexPeek(vexFile, verbose)
	if obsCode == 'Error':
		print 'Not processing file %s; it is not present or valid.' % vexFile
		exit(0)
	
	db = psycopg2.connect(databaseName)

	stations = stationTimes.keys()
	stations.sort()

	begin = 1.0e12
	end = -1.02e12
	for stn in stations:
		if stationTimes[stn][0] < begin:
			begin = stationTimes[stn][0]
		if stationTimes[stn][1] > end:
			end = stationTimes[stn][1]

	print 'Doing database query for all stations...'
	ccdata = getCableCalsFromDB(db, begin, end, verbose)

	for stn in stations:
		print 'Getting cable cal values for %s:' % stn

		outfile = lower('%s%s.%s.cablecal' % (obsCode, obsSeg, stn))
		if isfile(outfile) and not force:
			print 'Output file %s exists.  Won\'t overwrite.\n' % outfile
			continue

		out = open(outfile, 'w')
		out.write('# Cable cal data extracted from VLBA Monitor Database %s by %s ver. %s\n' % (split(databaseName)[0], program, version))
		out.write('# Ant MJD dur(days) C.Cal(psec)\n')

		for cc in ccdata:
			if cc[0] != stn:
				continue
			if cc[1] < stationTimes[stn][0] or cc[1] > stationTimes[stn][1]:
				continue
			# logic below modeled after old VLBA control code in module/getrtm.c
			# calc cf as shown on page 9 of VLBA Technical Report #5
			# 36.0/512.0 is the scale factor used by mon2xml
			# 3267.8 is a final scale factor used by mon2xml
			cin = cc[2]*512.0/36.0*3276.8
			cquad = cc[3]*512.0/36.0*3276.8

			# correct for mon2xml getting signedness confused.  These should be unsigned values.
			if cin < 0.0:
				cin += 16777216
			if cquad < 0.0:
				cquad += 16777216

			cf = cin / CMAX
			if (cin <  CHALF and cquad >= CHALFM) or (cin >= CHALF and cquad >= CHALFP):
				cf = 2.0 - cf
			picosec = cf * 500.0

			out.write('%s %14.8f 0.0 %10.6f  %14.10f %14.10f\n' % (upper(stn), cc[1], picosec, cc[2], cc[3]))

		out.close()

#---------

vexFile = ''
force = False
verbose = False

for a in argv[1:]:
	if a[0] == '-':
		if a in ['-f', '--force']:
			force = True
		elif a in ['-v', '--verbose']:
			verbose = True
		else:
			print 'Unknown option: %s' % a
			exit(0)
	else:
		if vexFile == '':
			vexFile = a
		else:
			print 'Unexpected parameter: %s' % a
			exit(0)

if vexFile == '':
	print 'No vex file provided.'
	exit(0)

processVexFile(vexFile, force, verbose)