python-snippets/SydneyTide-python.py

import datetime
import time
import math
from pylab import *

# Brett Miller  9th October 2015.
# This program reads Harmonics for Fort Denison as provided by BOM on the 9th October 2015
#
# NOT FOR DISTRIBUTION 
#   the predictions from this code have been cross checked against some daily predictions but 
#   not against all.
# USE AT YOUR OWN CAUTION.  NO GAURENTEE MADE THAT THESE PREDICTIONS ARE CORRECT.
#
# The x-tide method for calculating tides is:-
#  z(t,y,site)= MWL(site) + sum[c=1,N](nodefactor(c,y).amp(c,site).cos( speed(c).t + eqarg(c,y) - kappa(c,site) ))
#     where 
#	t is the time since the beginning of the year
#	y is the year
#	site is the site of interest
#	N is the number of consituents being included
#	nodefactor is an amplitude adjustment for each constituent for each year
#	eqarg is a phase adjustment for each site for each year
#	amp is the amplitude of each constituent for the site
#	kappa is the phase adjustment for each constituent for the site
#
# the harmonics included below in the python code were provided directly to WRL from BOM.


harm={}
speed=0
amp=1
lag=2

def readlinenocomm(F):
		l=F.readline()
		while l[0]=='#':
				l=F.readline()
		return l

#
H=open("HARMONIC","r")
l=readlinenocomm(H)
numconstits=int(l)

print "Number of Constituents = ",numconstits

# Read in the constituent speeds

constitspeed={}
i=0
while i<numconstits:
	l=readlinenocomm(H)
	ll=l.split()
	constitspeed[ll[0]]=float(ll[1])
	i=i+1
	
l=readlinenocomm(H)
sy=int(l)
l=readlinenocomm(H)
numyears=int(l)

# Read in the equalibrium arguements

equalibrium_arg={}
i=0
while i<numconstits:
	l=readlinenocomm(H)
	ll=l.split()
	constit=ll[0]
	v=[]
	j=0
	while j<=numyears/10:
		l=readlinenocomm(H)
		ll=l.split()
		for lll in ll:
			v.append(float(lll))
		j=j+1
	equalibrium_arg[constit]=v
	i=i+1
	
l=readlinenocomm(H)
if l[0:5]<>'*END*':
	print "Error, expected end after eq args.  Got:  ",l

# Read in the node factors
	
l=readlinenocomm(H)
numyears=int(l)

node_factor={}
i=0
while i<numconstits:
	l=readlinenocomm(H)
	ll=l.split()
	constit=ll[0]
	v=[]
	j=0
	while j<=numyears/10:
		l=readlinenocomm(H)
		ll=l.split()
		for lll in ll:
			v.append(float(lll))
		j=j+1
	node_factor[constit]=v
	i=i+1
	
l=readlinenocomm(H)
if l[0:5]<>'*END*':
	print "Error, expected end after node factors.  Got:  ",l

harm['sa']=[0.0410686,0.04,44.3373]
harm['ssa']=[0.0821373,0.0233,140.8642]
harm['Mm']=[0.5443747,0.0022,163.3925]
harm['MSF']=[1.0158958,0.0019,73.4003]
harm['MF']=[1.0980331,0.0025,203.8471]
harm['2Q1']=[12.8542862,0.0036,9.8664]
harm['SIGMA1']=[12.9271398,0.0043,16.6066]
harm['Q1']=[13.3986609,0.023,52.3723]
harm['RHO1']=[13.4715145,0.0043,55.9113]
harm['O1']=[13.9430356,0.0963,79.6431]
harm['MP1']=[14.0251729,0.001,45.9175]
harm['M1']=[14.4920521,0.0044,99.5325]
harm['CHI1']=[14.5695476,0.0013,102.4847]
harm['PI1']=[14.9178647,0.0028,125.9361]
harm['P1']=[14.9589314,0.0439,115.4037]
harm['S1']=[15,0.0033,129.8967]
harm['K1']=[15.0410686,0.148,119.6599]
harm['PSI1']=[15.0821353,0.002,11.0025]
harm['PHI1']=[15.1232059,0.0015,103.8567]
harm['THETA1']=[15.5125897,0.0017,137.3796]
harm['J1']=[15.5854433,0.0096,142.7051]
harm['SO1']=[16.0569644,0.0017,173.2765]
harm['OO1']=[16.1391017,0.0062,180.0375]
harm['OQ2']=[27.3416964,0.0006,75.4143]
harm['MNS2']=[27.4238337,0.0061,193.1893]
harm['2N2']=[27.8953548,0.0205,199.6233]
harm['MEU2']=[27.9682084,0.021,212.8684]
harm['N2']=[28.4397295,0.1127,223.5144]
harm['NEU2']=[28.5125831,0.0206,223.5377]
harm['OP2']=[28.9019669,0.0024,128.4082]
harm['M2']=[28.9841042,0.5022,237.061]
harm['MKS2']=[29.0662415,0.0014,313.3572]
harm['LAMDA2']=[29.4556253,0.0058,224.4572]
harm['L2']=[29.5284789,0.0135,243.1576]
harm['T2']=[29.9589333,0.0073,281.2802]
harm['S2']=[30,0.1247,260.7625]
harm['R2']=[30.0410667,0.0022,206.1633]
harm['K2']=[30.0821373,0.0373,250.3219]
harm['MSN2']=[30.5443747,0.0016,124.387]
harm['KJ2']=[30.626512,0.0023,62.9849]
harm['2SM2']=[31.0158958,0.0019,145.4123]
harm['MO3']=[42.9271398,0.0004,148.654]
harm['M3']=[43.4761563,0.0023,327.537]
harm['SO3']=[43.9430356,0.0002,307.6915]
harm['MK3']=[44.0251729,0.0005,248.2678]
harm['SK3']=[45.0410686,0.0016,113.2852]
harm['MN4']=[57.4238337,0.001,90.6269]
harm['M4']=[57.9682084,0.0031,114.2894]
harm['SN4']=[58.4397295,0.0003,103.5135]
harm['MS4']=[58.9841042,0.0014,183.7788]
harm['MK4']=[59.0662415,0.0004,196.3685]
harm['S4']=[60,0.0006,299.9646]
harm['SK4']=[60.0821373,0.0002,45.1703]
harm['2MN6']=[86.407938,0.0009,6.8059]
harm['M6']=[86.9523127,0.002,52.2368]
harm['MSN6']=[87.4238337,0.0007,64.4301]
harm['2MS6']=[87.9682084,0.0025,114.6867]
harm['2MK6']=[88.0503457,0.0007,131.372]
harm['2SM6']=[88.9841042,0.0003,157.8974]
harm['MSK6']=[89.0662415,0.0003,204.7339]
harm['2MN2S2']=[26.4079379,0.0003,1.3443]
harm['3M(SK)2']=[26.8701753,0.0005,3.6145]
harm['3M2S2']=[26.9523127,0.0008,78.7665]
harm['MNK2S2']=[27.505971,0.0001,225.556]
harm['SNK2']=[28.3575922,0.0003,218.9001]
harm['2SK2']=[29.9178627,0.0011,321.0943]
harm['MQ3']=[42.3827651,0.0004,43.8537]
harm['2MP3']=[43.009277,0.0001,134.0684]
harm['2MQ3']=[44.5695475,0.0001,329.7593]
harm['3MK4']=[56.8701754,0,39.731]
harm['3MS4']=[56.9523127,0.0001,97.9154]
harm['2MSK4']=[57.8860711,0.0002,283.9222]
harm['3MK5']=[71.9112441,0.0002,182.3575]
harm['M5']=[72.4602605,0,24.2117]
harm['3MO5']=[73.0092771,0.0005,264.7303]
harm['2(MN)S6']=[84.8476674,0,129.7631]
harm['3MNS6']=[85.3920422,0.0001,260.1927]
harm['4MK6']=[85.8542795,0.0001,332.2482]
harm['4MS6']=[85.9364168,0.0001,249.9087]
harm['2MSNK6']=[86.3258006,0,280.3885]
harm['2MNU6']=[86.4807915,0.0003,28.4884]
harm['3MSK6']=[86.8701754,0.0001,222.0912]
harm['4MN6']=[87.4966873,0.0003,292.3403]
harm['3MSN6']=[88.5125832,0.0006,338.3166]
harm['MKL6']=[88.5947204,0.0001,176.1921]
harm['2MN8']=[114.8476674,0.0001,219.9735]
harm['3MN8']=[115.3920422,0.0002,245.3801]
harm['M8']=[115.9364168,0.0002,284.0175]
harm['2MSN8']=[116.4079379,0.0002,203.8838]
harm['3MS8']=[116.9523127,0.0001,280.466]
harm['3MK8']=[117.03445,0.0001,286.3325]
harm['MSNK8']=[117.505971,0,219.0404]
harm['2MS8']=[117.9682084,0.0002,188.3974]
harm['2MSK8']=[118.0503457,0.0001,195.8549]
harm['4MS10']=[145.9364168,0.0002,293.7057]
harm['3M2S10']=[146.9523127,0.0002,307.5005]
harm['4MSN12']=[174.3761463,0.0001,110.9933]
harm['5MS12']=[174.920521,0,180.6601]
harm['4M2S12']=[175.9364168,0,108.791]
harm['SP3']=[44.9589314,0.0019,112.6443]
harm['S3']=[45,0.0012,23.4633]
harm['MA2']=[28.9430375,0.0039,117.2716]
harm['Ma2']=[29.0251709,0.0027,150.6304]
harm['MSV2']=[30.4715211,0,98.718]
harm['SKM2']=[31.0980331,0.0008,156.4092]
harm['2MNS4']=[56.407938,0.0001,26.7109]
harm['MV4']=[57.4966873,0.0003,94.1175]
harm['3MN4']=[58.5125832,0.0003,339.4208]
harm['2MSN4']=[59.5284789,0.0001,264.4342]
harm['NA2']=[28.3986628,0.0004,16.9727]
harm['MSO5']=[72.9271398,0.0001,247.7639]
harm['MSK5']=[74.0251729,0.0002,189.3859]

MSL=0.9952
		
# Now do the maths

start_time=datetime.datetime(2016,07,29,0,0,0)
end_time=datetime.datetime(2016,07,30,0,0,0)

time_step=datetime.timedelta(minutes=15)		

x=[]
y=[]
F=open("tide.out","w")
t=start_time
while (t<end_time):
	offset = t - datetime.datetime(t.year,1,1,0,0,0)
	yr=t.year-sy
	t1=(offset.days) * 24.0 + offset.seconds / 3600.0
	z=MSL
	for constit in harm.keys():
		if constit in node_factor:
			nf=node_factor[constit][yr]
			ea=equalibrium_arg[constit][yr]
		else:
			nf=1.0
			ea=0.0
		z=z+nf*harm[constit][amp]*math.cos((harm[constit][speed]*t1 + ea - harm[constit][lag]) * (2*3.1415927 / 360.0))

#	print t,z
	F.write("%s,%s\n" % (t,z))
	x.append(t)
	y.append(z)
	t=t+time_step

F.close
plt.figure()
plot(x,y,c='g',label='Sample')
plt.gcf().autofmt_xdate()
show()