NavList:
A Community Devoted to the Preservation and Practice of Celestial Navigation and Other Methods of Traditional Wayfinding
Some code for calculating tides
From: Frank Reed CT
Date: 2006 Jul 3, 17:37 -0500
For those of you who would like to experiment with tidal
calculations, here's some code I wrote a few years ago
for generating tide tables from a set of tidal harmonics. You
can get tidal harmonics for US locations from the NOS web
site. This code is written in Basic. I have marked the sections
that would need to be modified in a new implementation
with "##".
To see this code in action, visit www.fer1.com and click on Tides.
-FER
>>>>>
'some global variables:
dim pi, pio2, kk
dim TideH(25), rate(25), K(25), nodef(25), V0(25), u(25), harm(25)
dim lab(25), lin(25), hx(25), kx(25)
'TideH() is the array of tidal harmonic amplitiudes. K() is the array of
'the associated phases for each harmonic. rate() is the rate of each
'harmonic in degrees per hour.
'nodef() is the array of "node factors". As the location of the Moon's
'ascending node slowly precesses around the ecliptic in its 18 year cycle,
'the relative strength of the tidal harmonics changes.
'V0() is the array of initial phase angles at the epoch of the constants
'(the epoch date is Jan. 1, 2000).
'u() is the array of phases for the node factors.
'harm() is the array of standard names for the harmonics.
'## lab(),lin(),hx(),kx() are arrays for reading in the input data.
'pi and friends:
pi=3.14159
pio2 = pi / 2
kk = 180 / pi
'the total number of defined tidal harmonics:
MAXI = 17
'the standard names for the tidal harmonics:
harm(1) = "M2"
harm(2) = "N2"
harm(3) = "S2"
harm(4) = "K1"
harm(5) = "K2"
harm(6) = "O1"
harm(7) = "P1"
harm(8) = "M4"
harm(9) = "SA"
harm(10) = "SSA"
harm(11) = "L2"
harm(12) = "MU2"
harm(13) = "NU2"
harm(14) = "2N2"
harm(15) = "M6"
harm(16) = "Q1"
harm(17) = "T2"
'There are many more harmonics in the full theory, but this set will usually give times
'accurate to about five minutes and heights/currents accurate to a few percent.
'Weather creates larger variability than these limits.
'## insert code here to read in a tidal constant file. Here's an example:
'New London NOS
'zc: 5 'zc is the zone correction
'uni: feet 'uni are the units for the amplitudes.
'lat: 41.355 'lat/lon not strictly necessary
'lon: -72.087
'H0: 1.55 'H0 is the height of the tidal datum
'K1: 0.244 103.2
'K2: 0.060 280.3
'M2: 1.220 273.3 'Harmonics by name with amplitudes and phases.
'M4: 0.085 53.8 'These are derived from observations.
'N2: 0.285 252.4
'O1: 0.174 139.2
'P1: 0.074 115.0
'S2: 0.223 279.9
'SA: 0.181 135.0
'SSA: 0.059 68.9
'L2: 0.036 339.4
'MU2: 0.045 207.8
'NU2: 0.056 262.2
'2N2: 0.041 233.9
'M6: 0.039 131.3
'Q1: 0.050 126.4
'T2: 0.022 262.6
'## the lines below will need to be modified depending on the above
'code for reading in the tidal constants:
ZoneCorr = ...
units = ...
TideH(0) = ...
for j=1 to MAXI
TideH(j) = ...
K(j) = ...
next
'for tides relative to Mean Sea Level, ignore the tidal datum:
if MSLrequested then TideH(0)=0
'----hourly rates-----
'rate in degrees per hour of each tidal harmonic:
rate(1) = 28.9841042
rate(2) = 28.4397295
rate(3) = 30
rate(4) = 15.0410686
rate(5) = 30.0821373
rate(6) = 13.9430356
rate(7) = 14.9589314
rate(8) = 57.9682084
rate(9) = 0.0410686
rate(10) = 0.0821373
rate(11) = 29.5284789
rate(12) = 27.9682084
rate(13) = 28.5125831
rate(14) = 27.8953548
rate(15) = 86.9523127
rate(16) = 13.3986609
rate(17) = 29.9589333
'----initial angle----
'the initial phase in degrees of each tidal harmonic on Jan. 1, 2000:
V0(1) = 136.28
V0(2) = 7.83
V0(3) = 0
V0(4) = 9.36
V0(5) = 199.89
V0(6) = 127.07
V0(7) = 350.03
V0(8) = 272.57
V0(9) = 279.97
V0(10) = 199.95
V0(11) = 91.207
V0(12) = 272.74
V0(13) = 41.19
V0(14) = 239.38
V0(15) = 48.85
V0(16) = 358.62
V0(17) = 2.97
'## add code to input monthi, dayi, yeari or julian date directly...
jd = juldate(yeari, monthi, dayi)
'## add code to input starting time in hours, minutes...
timenow = (hournow + minnow/60)/24
ZoneCorr = ZoneCorr/24
TimeOffset = 0
'tide calculation follows...
'the number of hours from the epoch for the tidal constants:
dT = 24*(jd-2451545+timenow)
Node = 125.07 - 0.0529539*dT/24
Node = Node/kk
Perigee = 83.29 + 0.1114065*dT/24
Perigee = Perigee/kk
'a bunch of geometric factors for the varying orbital node:
inc = (23.455+5.145*COS(Node)+0.2625*(1-COS(2*Node)))/kk
xi = 12.10*sin(Node-0.225*sin(Node))
xinu2 = -2.14*sin(Node)
nup = -8.94*SIN(Node-0.153*SIN(Node))
xinua = 11.05*SIN(Node-0.248*SIN(Node))
nu = 13.14*SIN(Node-0.203*SIN(Node))
nupp2 = -17.77*sin(Node-0.0759*sin(Node))
PL2 = Perigee - xi/kk
coti2 = tan(pio2 - inc/2)
RL2 = kk*atn(sin(2*PL2)/(coti2*coti2/6-cos(2*PL2)))
tani2 = tan(inc/2)
tani2 = tani2*tani2
RAL2 = SQR(1-12*tani2*cos(2*PL2)+36*tani2*tani2)
'factors related to the Moon's orbital inclination w.r.t. the Earth's equator:
cosi2 = cos(inc/2)
sini2 = sin(2*inc)
sinis = sin(inc)^2
f78 = cosi2^4/0.9154
f227 = SQR(0.8965*sini2^2+0.6001*sini2*COS(nu/kk)+0.1006)
f75 = sin(inc)*cosi2^2/0.38
f235 = SQR(19.0444*sinis*sinis+2.7702*sinis*COS(2*nu/kk)+0.0981)
f215 = f78*RAL2
'the "node factors" for each harmonic.
nodef(1) = f78
nodef(2) = f78
nodef(3) = 1
nodef(4) = f227
nodef(5) = f235
nodef(6) = f75
nodef(7) = 1
nodef(8) = f78*f78
nodef(9) = 1
nodef(10) = 1
nodef(11) = f215
nodef(12) = f78
nodef(13) = f78
nodef(14) = f78
nodef(15) = nodef(8)*f78
nodef(16) = f75
nodef(17) = 1
u(1) = xinu2
u(2) = xinu2
u(3) = 0
u(4) = nup
u(5) = nupp2
u(6) = xinua
u(7) = 0
u(8) = 2*xinu2 'M4
u(9) = 0
u(10) = 0
u(11) = xinu2-RL2 'L2
u(12) = xinu2
u(13) = xinu2
u(14) = xinu2
u(15) = 3*xinu2 'M6
u(16) = xinua
u(17) = 0
'----the main calculation loop----
do
Tide = TideH(0)
dTide = 0
for i=1 to MAXI
Tphase = rate(i)*dT + V0(i)+u(i)-K(i)
Tphase = (Tphase - 360*int(Tphase/360))/kk
Tamp = TideH(i)*nodef(i)
Tide = Tide + Tamp*cos(Tphase)
dTide = dTide - rate(i)*Tamp*sin(Tphase)
next
'## Output result "Tide"
loop until done '##set done to some condition
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to NavList@fer3.com
To , send email to NavList-@fer3.com
-~----------~----~----~----~------~----~------~--~---
From: Frank Reed CT
Date: 2006 Jul 3, 17:37 -0500
For those of you who would like to experiment with tidal
calculations, here's some code I wrote a few years ago
for generating tide tables from a set of tidal harmonics. You
can get tidal harmonics for US locations from the NOS web
site. This code is written in Basic. I have marked the sections
that would need to be modified in a new implementation
with "##".
To see this code in action, visit www.fer1.com and click on Tides.
-FER
>>>>>
'some global variables:
dim pi, pio2, kk
dim TideH(25), rate(25), K(25), nodef(25), V0(25), u(25), harm(25)
dim lab(25), lin(25), hx(25), kx(25)
'TideH() is the array of tidal harmonic amplitiudes. K() is the array of
'the associated phases for each harmonic. rate() is the rate of each
'harmonic in degrees per hour.
'nodef() is the array of "node factors". As the location of the Moon's
'ascending node slowly precesses around the ecliptic in its 18 year cycle,
'the relative strength of the tidal harmonics changes.
'V0() is the array of initial phase angles at the epoch of the constants
'(the epoch date is Jan. 1, 2000).
'u() is the array of phases for the node factors.
'harm() is the array of standard names for the harmonics.
'## lab(),lin(),hx(),kx() are arrays for reading in the input data.
'pi and friends:
pi=3.14159
pio2 = pi / 2
kk = 180 / pi
'the total number of defined tidal harmonics:
MAXI = 17
'the standard names for the tidal harmonics:
harm(1) = "M2"
harm(2) = "N2"
harm(3) = "S2"
harm(4) = "K1"
harm(5) = "K2"
harm(6) = "O1"
harm(7) = "P1"
harm(8) = "M4"
harm(9) = "SA"
harm(10) = "SSA"
harm(11) = "L2"
harm(12) = "MU2"
harm(13) = "NU2"
harm(14) = "2N2"
harm(15) = "M6"
harm(16) = "Q1"
harm(17) = "T2"
'There are many more harmonics in the full theory, but this set will usually give times
'accurate to about five minutes and heights/currents accurate to a few percent.
'Weather creates larger variability than these limits.
'## insert code here to read in a tidal constant file. Here's an example:
'New London NOS
'zc: 5 'zc is the zone correction
'uni: feet 'uni are the units for the amplitudes.
'lat: 41.355 'lat/lon not strictly necessary
'lon: -72.087
'H0: 1.55 'H0 is the height of the tidal datum
'K1: 0.244 103.2
'K2: 0.060 280.3
'M2: 1.220 273.3 'Harmonics by name with amplitudes and phases.
'M4: 0.085 53.8 'These are derived from observations.
'N2: 0.285 252.4
'O1: 0.174 139.2
'P1: 0.074 115.0
'S2: 0.223 279.9
'SA: 0.181 135.0
'SSA: 0.059 68.9
'L2: 0.036 339.4
'MU2: 0.045 207.8
'NU2: 0.056 262.2
'2N2: 0.041 233.9
'M6: 0.039 131.3
'Q1: 0.050 126.4
'T2: 0.022 262.6
'## the lines below will need to be modified depending on the above
'code for reading in the tidal constants:
ZoneCorr = ...
units = ...
TideH(0) = ...
for j=1 to MAXI
TideH(j) = ...
K(j) = ...
next
'for tides relative to Mean Sea Level, ignore the tidal datum:
if MSLrequested then TideH(0)=0
'----hourly rates-----
'rate in degrees per hour of each tidal harmonic:
rate(1) = 28.9841042
rate(2) = 28.4397295
rate(3) = 30
rate(4) = 15.0410686
rate(5) = 30.0821373
rate(6) = 13.9430356
rate(7) = 14.9589314
rate(8) = 57.9682084
rate(9) = 0.0410686
rate(10) = 0.0821373
rate(11) = 29.5284789
rate(12) = 27.9682084
rate(13) = 28.5125831
rate(14) = 27.8953548
rate(15) = 86.9523127
rate(16) = 13.3986609
rate(17) = 29.9589333
'----initial angle----
'the initial phase in degrees of each tidal harmonic on Jan. 1, 2000:
V0(1) = 136.28
V0(2) = 7.83
V0(3) = 0
V0(4) = 9.36
V0(5) = 199.89
V0(6) = 127.07
V0(7) = 350.03
V0(8) = 272.57
V0(9) = 279.97
V0(10) = 199.95
V0(11) = 91.207
V0(12) = 272.74
V0(13) = 41.19
V0(14) = 239.38
V0(15) = 48.85
V0(16) = 358.62
V0(17) = 2.97
'## add code to input monthi, dayi, yeari or julian date directly...
jd = juldate(yeari, monthi, dayi)
'## add code to input starting time in hours, minutes...
timenow = (hournow + minnow/60)/24
ZoneCorr = ZoneCorr/24
TimeOffset = 0
'tide calculation follows...
'the number of hours from the epoch for the tidal constants:
dT = 24*(jd-2451545+timenow)
Node = 125.07 - 0.0529539*dT/24
Node = Node/kk
Perigee = 83.29 + 0.1114065*dT/24
Perigee = Perigee/kk
'a bunch of geometric factors for the varying orbital node:
inc = (23.455+5.145*COS(Node)+0.2625*(1-COS(2*Node)))/kk
xi = 12.10*sin(Node-0.225*sin(Node))
xinu2 = -2.14*sin(Node)
nup = -8.94*SIN(Node-0.153*SIN(Node))
xinua = 11.05*SIN(Node-0.248*SIN(Node))
nu = 13.14*SIN(Node-0.203*SIN(Node))
nupp2 = -17.77*sin(Node-0.0759*sin(Node))
PL2 = Perigee - xi/kk
coti2 = tan(pio2 - inc/2)
RL2 = kk*atn(sin(2*PL2)/(coti2*coti2/6-cos(2*PL2)))
tani2 = tan(inc/2)
tani2 = tani2*tani2
RAL2 = SQR(1-12*tani2*cos(2*PL2)+36*tani2*tani2)
'factors related to the Moon's orbital inclination w.r.t. the Earth's equator:
cosi2 = cos(inc/2)
sini2 = sin(2*inc)
sinis = sin(inc)^2
f78 = cosi2^4/0.9154
f227 = SQR(0.8965*sini2^2+0.6001*sini2*COS(nu/kk)+0.1006)
f75 = sin(inc)*cosi2^2/0.38
f235 = SQR(19.0444*sinis*sinis+2.7702*sinis*COS(2*nu/kk)+0.0981)
f215 = f78*RAL2
'the "node factors" for each harmonic.
nodef(1) = f78
nodef(2) = f78
nodef(3) = 1
nodef(4) = f227
nodef(5) = f235
nodef(6) = f75
nodef(7) = 1
nodef(8) = f78*f78
nodef(9) = 1
nodef(10) = 1
nodef(11) = f215
nodef(12) = f78
nodef(13) = f78
nodef(14) = f78
nodef(15) = nodef(8)*f78
nodef(16) = f75
nodef(17) = 1
u(1) = xinu2
u(2) = xinu2
u(3) = 0
u(4) = nup
u(5) = nupp2
u(6) = xinua
u(7) = 0
u(8) = 2*xinu2 'M4
u(9) = 0
u(10) = 0
u(11) = xinu2-RL2 'L2
u(12) = xinu2
u(13) = xinu2
u(14) = xinu2
u(15) = 3*xinu2 'M6
u(16) = xinua
u(17) = 0
'----the main calculation loop----
do
Tide = TideH(0)
dTide = 0
for i=1 to MAXI
Tphase = rate(i)*dT + V0(i)+u(i)-K(i)
Tphase = (Tphase - 360*int(Tphase/360))/kk
Tamp = TideH(i)*nodef(i)
Tide = Tide + Tamp*cos(Tphase)
dTide = dTide - rate(i)*Tamp*sin(Tphase)
next
'## Output result "Tide"
loop until done '##set done to some condition
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to NavList@fer3.com
To , send email to NavList-@fer3.com
-~----------~----~----~----~------~----~------~--~---
