#! /usr/bin/env python
# -*- coding: cp1252 -*-
'''
***************************************************************************************
                G E N E R A T E   U K O O A   P 1 / 9 0   F I L E
 
                          O B S - G E O M A R   T Y P E (this version is adapted
to Geomar scripts and formats of 2014 by nl)
                          

NOTE:
   Check any shot points that falls on midnight. If EIVA Survey Computer shot time 
   replacement by GPS Shot time means that date should be altered, this must be done
   manually.  

DESCRIPTION:   
   Merges information from two sources: a) The navigation files from the EIVA Survey
   computer, which holds records of Shot point (SP) number, UTM coordinates, date,
   time, depth and other parameters and b) accurate GPS based timing of each shot.
   
   Sequence of processing steps:
   
   A] The script applies an offset to EASTING & NORTHING coordinates in EIVA navigation
      file, in order to adjust for distance between seismic source and position reference
      point. Offset adjustment it based on treating shot points as vectors: First the 
      difference vector between current and previous SP is calculated. This vector is
      divided by its own length, obtaining a unit length vector. The unit length vector
      is then multiplied by the offset distance, yielding a correction vector that is
      subtracted from current SP. First SP is treated differently though, as there is no
      previous SP in that case; instead, next SP is used to obtain correction vector.  
   B] Collect corresponding GPS time stamp from the other file. Convert GPS time stamp to
      UTC. Check that timestamps in EIVA and GPS-file differs only by a user-defined
      amount.
   C] Generate UKOOA file.

INPUT FILES:
    Must be in same directory as this script. Name of input files stated in Constant
    section below.

OUTPUT FILE:
    Will be placed in same directory as this script.

UKOOA P1/90 HEADER, UTM ZONE, OFFSET DISTANCE:
    See 'Constant' section below.

EXTERNAL LIBRARY USED:
    UTM to Lat/Long conversion - download library from:
       http://www.pygps.org/#LatLongUTMconversion
    Either:
    a) On Linux: Unpack it, and run "python setup.py install" as root.
       In order to build Python modules you must first install 'python-devel'.
    b) Windows/Linux: Put "LatLongUTMconversion.py" in the same directory as this file.

PYTHON ON WINDOWS:
    This script was developed on Windows version of Python. It should run on Unix/Linux
    without any changes. Download Python for Windows from http://www.python.org
    N O T E :  In UKOOA file, line termination will be CR+LF on Win, and LF on Linux/Unix.
    
VERSION:

    Ver   Date         By         Description
    ----------------------------------------------------------------------------------
    0.6  08 Aug 2014  O.M.       LINE 367: Fixed error condition: If GPS shot time is less then
                                 LEAP_SECONDS after midnight, also the day part of the time
                                 variable will be decremented (note: we use Python's time object that also keeps
                                 info on year-month-day). Thus the comparison with nav data
                                 event time would fail (incorrectly). Fixed by checking 'day'
                                 part of variable, and adding 24 hours if necessary.
    0.5     Aug 2014  Nina       Adjusted UKOOA P1/90 output formatting so Geomar is content.
	0.4  23 Dec 2013  O.M.       2013 version.
    0.3   8 Jun 2010  O.M.       Can now use coordinates in all hemispheres.
    0.2   4 Jun 2009  O.M.       Changes:
                                   a) SP numbering made consecutive
                                   b) Much more plot QC output
                                   c) Fixed round off error if GPS fraction was > 0.995 s 
    0.1   1 Jun 2009  O. Meyer   Initial version
    ----------------------------------------------------------------------------------

                             Department of Earth Science
                                University of Bergen
                                     N O R W A Y

***************************************************************************************
'''

import LatLongUTMconversion
import os
from string  import split
from os.path import join
from datetime import datetime, timedelta
from time import strptime, strftime, gmtime
from cmath import *

CR_LF                   = '\x0d\x0a'


#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#++++
#++++           You can edit this section:


LINE_NAME               = '32'        # Max 6 characters, only number
UTM_Zone                = '33W'          # MUST also include letter!
WGS84                   = 23             # Datum constant, see LatLongUTMconversion file.

SOURCE_OFFSET           = 67.0            # Offset shot point positions by this amount.
GPS_LEAPSECONDS	        = 16              # Set to zero, if using UiB time stamp equipment (it outputs time in UTC, so no adjustment is needed)
                                          # 
TIME_DIFF_TOLERANCE     = 2.0             # In seconds. +/- tolerance on difference between Survey Computer 
                                          # shot event time stamp, and time stamp from Ashtech GG24.
SHOT_DISTANCE           = 200.0           # In meter
SHOT_DISTANCE_TOLERANCE = 3.0             # Shot interval distance tolerance, in meter

EivaLogFile             = 'OBS_2014-Storfjorden_event_data_line3-1-fixed.txt' 	# Name of EIVA Survey Computer log file, fixed (get rid of excess SPs)
GPS_TimestampFile       = 'p3-1-Storfjorden-OBS-2014-gps-timestamps-fixed.txt' 	# Name of GPS Time stamp file (fixed)
UkooaFile               = 'OBS-2014-P3-1-UKOOA-geomar.txt'  		    		# Name of output file
																	# NOTE: Plotting file also generated, with "_plot" extension
PlotFile                = UkooaFile.split('.')[0] + '_plot.txt'    	# You don't need to change this.


UKOOA_Header          = \
'''H0100 SURVEY AREA| STORFJORDEN, SVALBARD
H0101 SURVEY DETAILS| OBS-2014, PROFILE 3-1
H0102 VESSEL DETAILS| R.V. HAAKON MOSBY
H0200 SURVEY DATE| AUGUST 2014
H0201 FILE CREATED| %s
H0202 FILE VERSION| UKOOA P1/1990 - MODIFIED, SEE BELOW
H0300 CLIENT| OBS-2014
H0400 GEOPHYSICAL CONTRACTOR| UNIVERSITY OF BERGEN, NORWAY
H0500 POSITIONING CONTRACTOR| N/A
H0600 POSITIONING PROCESSING| O.M. UNIVERSITY OF BERGEN, DEPT OF EARTH SCIENCE
H0700 POSITIONING SYSTEM
H0700 SURVEY COMPUTER SOFTWARE| NAVIPAC, EIVA, DENMARK
H0700 GPS PRIMARY| SEAPATH 200
H0700 GPS SECONDARY| TRIMBLE 200D
H0700 GPS DIFF CORRECTION UNIT| FUGRO SEASTAR MOD 3510 LR
H0800 SHOTPOINT POSITION| CENTER OF SOURCE ARRAY
H1000 CLOCK TIME| UTC
H0900 OFFSET SHIP SYSTEM TO SP| %s M BEHIND
H1200 SPHEROID AS SURVEYED| A,1/F WGS-84 Spheroid 6378137.000 298.2572236
H1300 SPHEROID AS PLOTTED| A,1/F WGS-84 Spheroid 6378137.000 298.2572236
H1400 GEODETIC DATUM SURVEYED| WGS-84
H1500 GEODETIC DATUM PLOTTED| WGS-84
H1700 VERTICAL DATUM SL| ECHO SOUNDER
H1800 PROJECTION| 001 UTM NORTHERN HEMISPHERE
H1900 ZONE| %s NORTHERN HEMISPHERE
H2000 GRID UNIT|      1 INTERNATIONAL METER     1.00
H2001 HEIGHT UNIT|    1 INTERNATIONAL METER     1.00
H2002 ANGULAR UNITS| 1 DEGREES
H2600 ORIGINAL POSITIONS REFERENCED TO GRAVITYMETER, FRAME 38, 24 M FROM STERN
H2600 SOURCE DETAILS
H2600     OFFSET CENTRE OF SOURCE TO POS REF: 60 METER
H2600     OFFSET CENTRE OF SOURCE TO STERN:    60-24=36 METER
H2600     4 AIRGUNS, TOTAL VOLUME: 4800 CUIN
H2600     NO 1: BOLT 1500 1200 CUIN
H2600     NO 2: BOLT 1500 1200 CUIN
H2600     NO 3: BOLT 1500 1200 CUIN
H2600     NO 4: BOLT 1500 1200 CUIN
H2600     AIR PRESSURE NOMINAL: 138 BAR
H2600     SOURCE DEPTH NOMINAL: 6 METER
H2600 MAGNETOMETER DETAILS
H2600     MFR: GEOMETRICS, MOD: G-882, SN: 882150
H2600     TOW FISH CABLE LENGTH DEPLOYED: 195 METER
H2600     OFFSET TO SOURCE: 195-36 = 159 METER
H2600     OFFSET TO NAVIGATION SYSTEM REFERENCE: 195+24 = 219 METER
H2600 MINI-STREAMER DETAILS
H2600     MFR: GEOMETRICS, MOD: MICROEEL
H2600     CHANNELS: 8, GL: 6.25 M, 50 M ACTIVE, SN: ME-0008
H2600     STREAMER DEPTH NOMINAL: NO DEPTH CONTROL
H2600     OFFSET MIDPOINT FIRST CH TO POS REF:  X M 
H2600     OFFSET MIDPOINT FIRST CH TO SOURCE:  X-Y=Z METER
H2600 GRAVITYMETER DETAILS
H2600     MFR: L&R (SENSOR AND FRAME) ZLS (ELECTRONICS), SN: S-99
H2600     SOFTWARE VERSION: 1.98 (FOGPACK GYROS)
H2600
H2600 FILE FORMAT MODIFIED TO INCLUDE SHOTPOINT TIME IN COLUMNS [7..19]
H2600
H2600 00011111111112222222222333333333344444444445555555555666666666677777777778
H2600 78901234567890123456789012345678901234567890123456789012345678901234567890
H2600 <--SP TIME--><-SP-><---LAT--><--LONG---><--EAST-><-NORTH-><-DEP><D><TIME>-''' % (strftime("%d %b %Y %H:%M:%S UTC", gmtime()).upper(),
                                                                                         SOURCE_OFFSET,
                                                                                         UTM_Zone
                                                                                         )


#++++
#++++           End of user edit section
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++





#----------------------------------------------------------------------------------------------------------
#  -----     S U B R O U T I N E S
#----------------------------------------------------------------------------------------------------------


#--------------------------------------------------------------------------------------------------
def Degrees2DMS(decimal_degree, lat_lon):     #----- Convert from decimal degrees to dddmmss.ss.
#--------------------------------------------------------------------------------------------------

    hemisphere = ""
    lat_lon    = lat_lon.upper()
    
    
    if lat_lon == "LAT":
        if decimal_degree > 90.0:
            print "Error: LAT value > 90.0 deg:", decimal_degree
            os.exit()
        hemisphere = "N"             # Assume Northern hemisphere
        if decimal_degree < 0.0:
            hemisphere = "S"                            # Wrong assumption ...
            decimal_degree = abs(decimal_degree)        # Get rid of '-'
    
    
    elif lat_lon == "LON":
        if decimal_degree > 180.0:
            decimal_degree = 360.0 - decimal_degree     #----- Transformation: [0..360> -> <-180..180]
        
        hemisphere = "E"                                # Assume Eastern hemisphere
        if  decimal_degree < 0.0:
            hemisphere = "W"                            # Wrong assumption ...
            decimal_degree = abs(decimal_degree)        # Get rid of '-'
            
            
    if not hemisphere:
        print "Error: Lat/Lon must be specified by either 'LAT' or 'LON' (small/cap letters accepted)." 
        os.exit()

    
    decideg = str(decimal_degree)
        
    s = decideg.split('.')
    degree_main = s[0]
    degree_fraction = float(decideg) - float(degree_main)
    minute = str(degree_fraction * 60)
    minute_main = minute.split('.')[0]
    minute_fraction = '0.' + minute.split('.')[1]
    sec = float(minute_fraction) * 60
    second = "%5.2f" % sec

    return("%s%2d%s%s" % (degree_main, int(minute_main), second, hemisphere))



#---------------------------------------------------------------------
#-----            S t a r t   o f   p r o g r a m
#---------------------------------------------------------------------


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>  1 s t   s e c t i o n  
#>>
#>>  Loop through all lines in EIVA shot point file, and extract information.
#>>  Convert UTM to Lat, Lon. Store in dictionary for later use.
#>>  Get corresponding entry in GPS Time stamp file, convert to UTC, compare to 
#>>  shot time in EIVA Survey Computer log file. Exit if errors are found.
#>>
#>>  Example line from Eiva shot point file:
#>>  000057 2010.06.02 02:51:41.914 0334654.77 07118140.71 00.90 038.1 00000.15 005.14 000.000 "H Mosby f " 3002.22 9449.08 52045.27 "" 01
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


f           = open(EivaLogFile)
GPS_file    = open(GPS_TimestampFile)
GPS_LineNumber  = 0

Eiva_Dictionary          = {}               # Will store all info extracted here. SP number used as key.
GPS_Timestamp_Dictionary = {}

LeapSeconds         = timedelta(seconds=GPS_LEAPSECONDS)
TimeDiffTolerance   = timedelta(seconds=TIME_DIFF_TOLERANCE)

SeekFirstSP = True                  		# Ensure SP numbering is consecutive

for line in f:
    if line[0] == '0':              		# Assume SP is first, with leading zeros. Gets rid of comments, header this way.
        s = line.split()            		# Split line into array    
        if SeekFirstSP:
            Shotpoint = int(s[0])
            SeekFirstSP = False
        else:
            Shotpoint += 1
        Date      = s[1]
        NavTime   = s[2]
        Easting   = float(s[3])
        Northing  = float(s[4])
        Dol       = float(s[8])
        Depth     = float(s[14])

        (Lat, Lon) = LatLongUTMconversion.UTMtoLL(WGS84, Northing, Easting, UTM_Zone)       # Convert EASTING & NORTHING to lat,lon
        
        #print Shotpoint, Easting, Northing, Lat, Lon
        #z = raw_input("Proceed? ")

        #---- Use date & time info to declare Python time object.
        #---- Using these time objects makes it easy to perform time adjustments.

        Year     = Date.split('.')[0]
        Month    = Date.split('.')[1]
        Day      = Date.split('.')[2]
        Hour     = NavTime.split(':')[0]
        Minute   = NavTime.split(':')[1]
        Second   = NavTime.split(':')[2].split('.')[0]
        Microsec = NavTime.split(':')[2].split('.')[1]

        EivaTimeStamp = datetime(int(Year),
                                 int(Month),
                                 int(Day),
                                 int(Hour),
                                 int(Minute),
                                 int(Second),
                                 (int(Microsec) * 1000)
                                 )

        
        #print Shot point, EivaTimeStamp, Correction, NewTimeStamp                            # Uncomment if you would like to check how it works
        
        #--- Store info in dictionary (with SP as key):    
        Eiva_Dictionary[Shotpoint]  = (Easting, Northing, EivaTimeStamp, Depth, Lat, Lon, Dol)
        
        #print Shot point, Date, Time, EivaTimeStamp, Easting, Northing, Lat, Lon, Depth      # Uncomment if you like watching scrolling text ...
        
        
        # -------------------------------------------------------------------------------
        #---- Extract GPS time and subtract number of leap seconds (16, as of Sept 2013)
        #---  GPS time stamp example
        #       $PASHR,TTT,4,00:59:48.7383548*0C      
        #       $PASHR,TTT,4,01:01:13.9997177*07
		
		#       $PASHR,TTT,1,03:48:52.0042740*02      # 2013, new Ashtech unit
		#       $PASHR,TTT,1,03:50:19.7022105*02
        
        
        # Skip GPS time stamp lines starting with # (i.e. comments)
        while True:
            GPS_line = GPS_file.readline()
            GPS_LineNumber += 1
            if GPS_line[0] != "#":
                break

        try:
            t = GPS_line.split(',')
            timestamp = t[3][:-5]     					# Cut five last characters in time stamp (CR + checksum + last digit)
            print "GPS time:", timestamp,				# ENSURE that resulting fraction has 6 digits, i.e. represents microseconds
														# If we round microseconds we must take care of seconds etc. also
            hh = timestamp.split(':')[0]
            mm = timestamp.split(':')[1]
            ss = timestamp.split(':')[2]
            print '>',

            sec = ss.split('.')[0]                     # Split seconds in whole part and ...
            us  = ss.split('.')[1]                     #    fraction (microseconds)
            print '(sec=%s, us=%s)  subtract: %d leap seconds' % (sec, us, GPS_LEAPSECONDS),
            if len(us) != 6:
                print
                print '---->', GPS_line
                print 'Error: microsecond value MUST have 6 digits.'
                f.close()
                GPS_file.close()
                exit(1)
            time_GPS = datetime(2013,12,23,int(hh),int(mm),int(sec),int(us))        # Load DateTime Object. Use any values for year, month, date
            print '=>',
            
            time_UTC = time_GPS - LeapSeconds                                       # Subtract leap seconds to get UTC time

            if (time_GPS.day !=  time_UTC.day):                 # If we rolled back to previous day, add 24 hours
                time_UTC += timedelta(hours=24)                 # otherwise the comparison with nav time will FAIL
            
        except:
            print
            print '---->', GPS_line
            print 'Error parsing this line in GPS timestamp file....'
            print "Line number in GPS time stamps file:", GPS_LineNumber
            f.close()
            GPS_file.close()
            exit(1)
        
        ###if len(GPS_line) != 33:
        ###    print '---->', GPS_line
        ###    print 'Error, length of this line is not 33 characters ...'
        ###    f.close()
        ###    GPS_file.close()
        ###    exit(1)
      

        #--- Check if round-off to milliseconds affects seconds (and possible higher order also):

        if time_UTC.microsecond >= 999500:
            time_UTC += timedelta(seconds=1)                              # Add one second - will take care of higher order also!
            time_UTC -= timedelta(microseconds=time_UTC.microsecond)      # Get rid of second fraction
                                                                           # From now on, UTC time can be rounded to milliseconds without affecting higher parts
        
        #--- Now format UTC time string, also round off microseconds to milliseconds

        us   = '%1.3f' % (time_UTC.microsecond/1000000.0)
        us   = us.split('.')[1]
##nl        time = '%02d%02d%02d.%s000' % (time_UTC.hour, time_UTC.minute, time_UTC.second, us)
        time = '%02d%02d%02d.%s' % (time_UTC.hour, time_UTC.minute, time_UTC.second, us)
        print "UTC:", time

        time_QC = "%s %s" % (t[3][:-4], time)   ### Use this for Ashtech GG24 data
        ###time_QC = "%s %s" % (t[3], time)           ### Use this for UiB Time stamp unit
        
        GPS_Timestamp_Dictionary[Shotpoint] = time_QC        # Store in dictionary (key = SP): String containing UTC time, GPS time stamp

 



        
        #-------------------------------------------------------------------------------------------------
        #---  Check: Is there match between shot point times, one from Survey Computer log file, the other
        #---         from GPS time stamp (and leap seconds have been subtracted, so they should
        #            match pretty good). Tolerance is set in constant section in top of this file.
        #            Time in navigation data log file assumed to have this format: 17:10:41.306

        # This section is also a bit messy, clean up, please .... (but it works)

        hh_navdata = NavTime.split(':')[0]
        mm_navdata = NavTime.split(':')[1]
        ss_navdata = NavTime.split(':')[2].split('.')[0]        # Discard fractions of seconds ..
        
        # Get microseconds by subtracting whole part of second from total second value, multiply
        # the resulting fraction by 10e6, and cast to int
        
        us_navdata = int((float(NavTime.split(':')[2]) - float(ss_navdata)) * 1000000)
        # Use any values for year, month, date ...
        t_navdata = datetime(2013,12,23,int(hh_navdata),int(mm_navdata),int(ss_navdata),int(us_navdata))

        UpperLimit = t_navdata + TimeDiffTolerance
        LowerLimit = t_navdata - TimeDiffTolerance
        if (time_UTC > UpperLimit) or (time_UTC < LowerLimit):		# tc is GPS time (leap second adjusted)
            print
            print '---> Nav data time:', t_navdata
            print '---> Upper limit:', UpperLimit
            print '---> Lower limit:', LowerLimit
            print '---> Adjusted GPS time:', time_UTC
            print '--->', line[:35]
            print '---> GPS time stamp (adjusted for leap seconds):', time
            print '---> Mismatch Survey Computer time stamp and GPS time stamp of shot event'
            print "Line number in GPS time stamps file:", GPS_LineNumber
            f.close()
            GPS_file.close()
            exit(1)
        

f.close()
GPS_file.close()

##exit(0)





#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>           2 n d   s e c t i o n         
#>>
#>>  Offset adjustment to seismic source
#>>  Using vector calculus and let complex numbers represent vectors.
#>>  So each Shot point is a vector -> 
#>>      EASTING=real part, NORTHING=imaginary part of complex number
#>>
#>>  Complex numbers prefixed by "c" (like cThisPosition)
#>>
#>> Also writes plot file (to compare old/new positions). 
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


f = open(PlotFile, "w")
f.write("#SP, oldEasting, oldNorthing, oldLat, oldLon, newEasting, newNorthing, newLat, newLon, Offset, deltaPos, deltaTime, Speed, Dol, time_GPS, time_UTC\n")

First_SP = True

SP_keys = Eiva_Dictionary.keys()
SP_keys.sort()

PositionDictionary = {}         #   Will use SP as key in this dictionary
DeltaTime = timedelta()

for SP in SP_keys:              #<<<--- Loop through every shot point
    (ThisEasting, ThisNorthing, ThisNewTimeStamp, ThisDepth, ThisLat, ThisLon, ThisDol) = Eiva_Dictionary[SP]       # Retrieve info from dictionary, SP = key
    cThisPosition            = ThisEasting + 1j*ThisNorthing
    
    if not First_SP:
        (PrevEasting, PrevNorthing, PrevNewTimeStamp, PrevDepth, PrevLat, PrevLon, PrevDol) = Eiva_Dictionary[SP-1] 
        cPrevPosition        = PrevEasting + 1j*PrevNorthing
        cDeltaPosition       = cThisPosition - cPrevPosition
        DeltaPos             = abs(cDeltaPosition)
        if (abs(DeltaPos - SHOT_DISTANCE) > SHOT_DISTANCE_TOLERANCE):
            print 'Warning SP=%d: Delta Pos %.1f m (%0.1f +/-%.1f m)' % (SP, DeltaPos,SHOT_DISTANCE, SHOT_DISTANCE_TOLERANCE)
        DeltaTime = ThisNewTimeStamp - PrevNewTimeStamp
        DeltaTime_f = (float(DeltaTime.seconds) + (float(DeltaTime.microseconds) / 1000000) )
        Speed = DeltaPos / DeltaTime_f     # in m/s
        Speed = Speed * 3600/1852          # in knots
        ### print DeltaTime, "%.1f" % DeltaTime_f, Speed
    else:   # This was the first SP, must treat it differently
        First_SP = False
        (NextEasting, NextNorthing, NextNewTimeStamp, NextDepth, NextLat, NextLon, NextDol) = Eiva_Dictionary[SP+1]
        cNextPosition        = NextEasting + 1j*NextNorthing
        cDeltaPosition       = cNextPosition - cThisPosition
        DeltaPos             = 0.0
        Speed                = 0.0
        
    cOffsetVector            = SOURCE_OFFSET * cDeltaPosition/abs(cDeltaPosition)   # Offset in m multiplied by unit vector with proper direction
    cOffsetCorrectedPosition = cThisPosition - cOffsetVector

    #--- Transform new position to geographical coordinates.
    e = cOffsetCorrectedPosition.real
    n = cOffsetCorrectedPosition.imag
    (NewLat, NewLon) = LatLongUTMconversion.UTMtoLL(WGS84, n, e, UTM_Zone)
    
    PositionDictionary[SP]   = (cOffsetCorrectedPosition, NewLat, NewLon)
    #print SP, cThisPosition, cOffsetVector, abs(cOffsetVector), cOffsetCorrectedPosition       # Uncomment to see complex numbers in action
    
    time_QC = GPS_Timestamp_Dictionary[SP]
    
    # Write original and offset adjusted positions to plot file
    f.write("%s %.1f %.1f %.10f %.10f %.1f %.1f %.10f %.10f %.1f %.1f %.1f %.1f %.1f %s\n" % ( SP,
                                                                               cThisPosition.real,
                                                                               cThisPosition.imag,
                                                                               ThisLat,
                                                                               ThisLon,
                                                                               cOffsetCorrectedPosition.real,
                                                                               cOffsetCorrectedPosition.imag,
                                                                               NewLat,
                                                                               NewLon,
                                                                               abs(cOffsetVector),
                                                                               DeltaPos,
                                                                               (float(DeltaTime.seconds) + (float(DeltaTime.microseconds) / 1000000)),
                                                                               Speed,
                                                                               ThisDol,
                                                                               time_QC
                                                                               ))
    
f.close()
print
print 'Plot file: ' + PlotFile



##exit(0)




#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>          3 r d   s e c t i o n      
#>>
#>>  Write UKOOA P1/90 lines to file
#>>  Example start of data section:
#>>  
#>>     H2600 00011111111112222222222333333333344444444445555555555666666666677777777778
#>>     H2600 78901234567890123456789012345678901234567890123456789012345678901234567890
#>>     H2600 <--SP TIME--><-SP-><---LAT--><--LONG---><--EAST-><-NORTH-><-DEP><D><TIME>-
#>>     STSD-1132634.376000    47625937.91N  657 8.45E 396272.46986557.1  86.8148132634 
#>>     STSD-1132806.349000    48625942.26N  65657.96E 396129.06986696.6  89.9148132806 
#>>     STSD-1132926.968000    49625946.54N  65647.31E 395983.56986833.7  89.3148132926 
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


f = open(UkooaFile, "wb")

HeaderLines = UKOOA_Header.split('\n')   # --- First get UKOOA P1/90 header

Buffer = ''
NewLine = ''

for line in HeaderLines:
  s = line.strip()
  if '|' in s:
    t = s.split('|')
    try:
        NewLine = '%-32s%-48s\r\n' % ( t[0], t[1].strip() )
    except:
        NewLine = '%-32s%-48s\r\n' % ( t[0], '')
  else:
    NewLine = '%-80s\r\n' % s
    
  if len(NewLine) <> 82:
    print "Error in UKOOA header: Line not 82 characters long (incl <CR>+<LF>)"
    print "   ", NewLine
    exit(1)
  Buffer += NewLine
  
  

for SP in SP_keys:      # --- Loop over sorted list of shot point numbers

    #----  Fetch information from two directories, using shot point number as key
    
    (OldEasting, OldNorthing, EivaTimeStamp, Depth, OldLat, OldLon, Dol) = Eiva_Dictionary[SP]
    (cOffsetCorrectedPosition, Lat, Lon)                                 = PositionDictionary[SP]
    time_QC  = GPS_Timestamp_Dictionary[SP]
    Easting  = cOffsetCorrectedPosition.real
    Northing = cOffsetCorrectedPosition.imag
    
    time   = time_QC.split()[1]                      # Fetch (from dictionary) UTC time value

    ukooa  = 'S'                                      # First char of record = S, meaning "Centre of source"
##nl    ukooa += "%-5s"  % LINE_NAME                      # -5s means string, width=5, left justify
    ukooa += "%-6s"  % LINE_NAME                      # -5s means string, width=6, left justify
    ukooa +=  time                                    # Here we put OBS time stamp (not according to standard P1/90)
##nl    ukooa += "%6d"   % SP                             # Event = SP; width = 6, right justified
    ukooa += "%8d"   % SP                             # Event = SP; width = 8, right justified ##nl
#    ukooa += ' '                                      # Spare field, width = 1  ##nl
    ukooa += "%10s"  % Degrees2DMS(Lat, "LAT")        # Latitude
    ukooa += "%11s"  % Degrees2DMS(Lon, "LON")        # Longitude
    ukooa += "%9.1f" % Easting                        # Easting, 9 characters in total, one decimal point
    ukooa += "%9.1f" % Northing                       # Northing, --"--
    ukooa += "%6.1f" % Depth                          # Depth
    ukooa += "%3d"   % EivaTimeStamp.timetuple()[7]   # Day of year (using EIVA time data - could mean trouble at midnight - check manually)
##    ukooa += ' '                                      # Spare field, width = 1  ##nl
    ukooa +=  time.split('.')[0]                      # Time, without fraction (no rounding!)
##nl    ukooa += ' '                                      # Spare field, width = 1
    ukooa +=     '\r\n'                                     # This will generate CR+LF on Windows ....
              
    Buffer += ukooa
    if len(ukooa) > 82:
        print 'ERROR: UKOOA data line exceeding 82 characters'
        print ukooa,
        exit(1)
    #print ukooa,

##nl Buffer += 'EOF\r\n'

f.write(Buffer)

f.close()
print 'UKOOA P1/90 file: ' + UkooaFile
print "CHECK - did GPS timestamp correction cause roll-over at midnight?"
print "REMOVE the last empty line in the file"