Stormy Weather SoftWareTIDAL PREDICTION PROGRAMME NOTESThese notes are purely technical, will probably interest very few people, prove uninteligible to non  mathematicians and therefore might never be read. They are included for completeness of our Tidal Prediction SoftWare Some Manufacturers' Trade Marks are included (Microsoft, OS/2, Apple etc) and should be recognized as such. HISTORY.Version 1 was a virtually unusable set of math routines, based on various source documents, mainly by Doodson. User friendliness was nonexistant, but the results were usable, demonstrating the validity of "new" algorithms that we had developed specifically for computerised solutions after being prompted by Mike Harris (Pangolin Software) looking for simple but fast results from programmable calculators and very small computers. The new algorithms were adapted for use on PC type computers running a minimum of DOS 3.1, 80286 CPU and .5 Mb hard disk space available. Version 2 developed ease of use, after the decision was made to use sets of constituents containing a maximum of ten harmonics with the possibility of inferring six more. Seven nodal corrections were used. Version 2.1 followed a number of most useful contacts with the Canadian Hydrographic Service, and Version 2.2 (the latest fully usable version as of November 1995) responded to comments made by user/mariners to whom working copies had been made available. Accuracy, compared to published tide tables, is very acceptable; for semidiurnal and diurnal tides where little shoal water effect is felt errors rarely exceed 2% in height and 5 minutes in time. However mixed tides in estuaries can produce errors up to 8% in height and 15 minutes. Average errors are obviously lesser, and the majority of time errors occur at stands. Speed was improved to the point that a complete month of High/Low predictions was tabulated, ready to print, and produced graphically on screen within 15 seconds on a 286, and nearly instantaneously on a 486 DX (the video driver being a major "slow" point!) It has become evident that further development in accuracy is not possible without further harmonic input; also 32 bit CPU technology is developing so fast (OS/2, Windows 95, etc) that programming should now be orientated towards capability of compiling for DOS as well as 16 and 32 bit MS Windows. Minor variations for OS/2 are easily incorporated at compilation time. Current programming in QuickBasic 4.5, Basic 7.0 and VisualBasic 3.0 is adaptable to these ends. Investigations into a Macintosh/Apple version continue. If a C++ version becomes necessary, it will demand a complete rewrite. Fortran has been considered for the math routines, but the interface with other languages poses considerable difficulty. VERSION 3.x(References in this section are, alphabetically
Major changes include:Data File Structure:To include IHO reference system and 37 constituents (see below) in binary format. A fixed length record of 184 Bytes is required for each station. This is just about equal to v2.x ASCII records containing only 10 constituents. This is operational, together with utilities to convert other formats (Canada, US, UK, France) automatically. Astronomical Data:A complete new utility has been written based J2000. Its accuracy, after comparison with the Astronomical Almanach, is far superior to tidal requirements. Using a fixed reference date (1 Jan 1980 for v2.1 or 1 Jan 2000 for v2.2) with daily increments can, due to variations in orbital phenomena, plus Delta T to translate from TBT to UT, lead to errors well in excess of 1°. It seems a shame to use harmonic data correct to 0.1° under these circustances, even if the original analysis contained the same degree of uncertainty. Therefore it is proposed that daily increments from 0 hours, 16th day of each month be used. An added bonus is that single precision (7 sig figs) can be used in lieu of double (15 sig figs). Either a data file containing say 1995 to ? could be used, or the Tide programme could access the astro utility directly. Constituents:A total of 37 constituents can be used if available. See Appendix 1. The choice of 4 third diurnals and 3 sixth diurnals represents a major departure from v2.x. MF reviewed the list and suggested the inclusion of M3. Eighth diurnals and further have been neglected, therefore double highs (ex Southampton) may still not be quite accurate. 16 satellites (nodal corrections) have been retained at this point, pertaining to the main constituents. Selection was based on the Cartwright coefficient being greater than 0.1. Third order harmonics have been avoided where possible (quarter phase shifts and corrections for latitude), negative amplitudes can be corrected by a half phase shift. All Doodson numbers (s, h, p, N and p') are corrected by adding {Species x ( s +h )} (Species being 0 for long term, 1 for diurnals, .... 6 for sixth diurnals); this facilitates conversion from Tau (lunar day) to UTC in later calculations. Formulae:(1) ASTRO Rebased 1/1/2000 (Brackets 1980) where NumDay is the number of days (+ or ) from reference date and the function Modulus 360 maintains orbital data in range 0° to 360° s = FNMod360#(211.72826# + 13.17639673# * NumDay&) 'Long moon (78.15194) h = FNMod360#(279.86807# + .98564734# * NumDay&) 'Long Sun (279.72234) p = FNMod360#(83.29754# + .11140408# * NumDay&) 'Perigee (349.49439) N = FNMod360#(234.98196# + .05295392# * NumDay&) 'RA Moon (208.15463) p' = FNMod360#(282.93732# + .000047069# * NumDay&) 'Perihelion (282.59344) the result being for Time = 0 hours on NumDay (T0), partial times if required are found by linear interpolation; T24 = T0 on NumDay + 1 (2) SUM (Sigma) OF CONSTITUENTS by species where A = Amplitude, g = total phase lag at T0, repeated for T24. a) for each species (0, 1 .....6), starting SigmaX = 0, using major constituents and satellites SigmaCosine = SigmaCosine + A * Cosine(g) SigmaSine = SigmaSine + A * Sine(g) b) reduce to average harmonic for species where Ex is amplitude, Phi is phase lag Ex = SquareRoot ( SigmaSineý + SigmaCosineý ) [If Ex = 0 OR Ex = SigmaCosine then Ex = Ex + 10 ^ 10 to avoid division by 0] Phi = ArcCosine(SigmaCosine / Ex) [* 1 if SigmaSine is negative] We now have a series of "Average" amplitudes aand phases for each species at Time = 0 and Time = 24 for the day: Ex0 and Ex24, Phi0 and Phi24 BY SPECIES 3) HEIGHT OF TIDE at any Time (TX) Let suffixes 0, X and 24 represent Time in hours on any date, and maintaining 180° < (Phi24  Phi0) < 180°, then
4) EXTREMA found by using derivative for slope = 0 of type Sine/Cosine. The value of the derivative multiplied by K = 12.42 or 24.84 / 2 * Pi can be used to approach the absolute value of the Timestep (see 5 below) Using Formula Height in 3) as FHtCos, and rewriting with Sine in lieu of Cosine as ............... FHtSin d = K * { FHtSin(Species1) + 2*FHtSin(Species2) ....... 6*FHtSin(Species6)} _____________________________________________________________________ { FHtCos(Species1) + 4*FHtCos(Species2) ....... 36*FHtCos(Species6)} Note that integers have been used, but could be refined using average angular speed. 5) TIME STEPS for extrema are currently 12.42, 6.21 and 1.03 hours respectively for diurnal, semidiurnal and mixed tides, and works well as a "DeltaTime", to or from which the derivative above is added or subtracted as refinement in a recursive loop. The number of recursions averages 3.7 for a final "DeltaTime" < 1 minute. This will be modified to respect GG, and thus to ensure compatibility with CHS tables. 6) Further precision can be obtained (at the expense of computer time) by using times of extrema from 4), and returning for new values in steps 1) through 4) of each extremum. This does not appear to be necessary for pure diurnal or semidiurnal tides, but does show an improvement for mixed tides, particularly those that within a calendar month produce a mix of two or four extrema per day. However further investigation is necesary, above all in assessing the influence of shoal water effects. 7) Stands and double tides. These are currently assessed on the basis of height alone, and not recorded as extrema although they appear in the graphical daily representation (heights calculated for every 0.25 hours). The programme has now been modified so that 6 extrema per day can be recorded and will be refined, again along the lines suggested by GG, so that time as well as height (proposed as 0.01 * Z0) will be considered. Preliminary trials using Trois Rivières looks promising. Appendix 1: ConstituentsEach file record contains the following:
Doodson Numbers given as programmed (CHS, MF equivalents in brackets difference being {Species * (s  h)} Angular speed in degrees/solar hour (°/hr) to 7 sig figs Z0  Mean level  Constant Long Period, total = 5, Species 0:
Diurnals, total = 12, Species 1:
Repeats: O1=MK1 K1=MO1=SP1 SIG1=SIGMA1 Q1=NK1 J1=MQ1 M1=NO1 P1=SK1 MP1=TAU1 KQ1=UPS1 M(SK)2=H1 ? Series1 ? M(SK)2=H1 Semi Diurnals, total = 9, Species 2:N2 3s +2h +p s +p 28.43973Mu2 4s +4h 2s +2h 27.96821 L2 s +2h p 180° s p 180° 29.52848 SH s +2h p Amplit * 1 T2 h +p' 2s 3h +p' 29.95893 2N2 4s +2h +2p 2s +2p 27.89535 Nu2 3s +4h p s +2h p 28.51258 K2 2h 2s 30.08214 M2 2s +2h none 28.9841 S2 none 2s  2h 30.0 Not included : Lda2 OQ2 R2 H2 Gam2 Eps2 P2 MKS2 MSN2 Eta2 2SM2 OP2 2SM2 Repeats: M2=KO2 S2=KP2 N2=KQ2 L2=2MN2=L2 2MN2(Fr) 2N2=2MS2 KJ2=ETA2 MNS2=UPS2 M(KS)2=H2 ? Series2 ? 2MN2S2=ST36 M(SK)2=H1 Third Diurnals, total = 4, Species 3:M3 3s +3h 180° 180° 43.47616 AD M1+M2 = same +90° NOTE 2MO3 4s +3h 90° s 90° 42.92714 SO3 2s +h 90° s 2h 90° 43.94304 MK3 2s +3h 270° s 270° 44.02517 AD M2 + K1 = same +90° Not included: SK3 Repeats: MO3=2MK3 MQ3=NO3 Quarter Diurnals, total = 4, Species 4:M4 4s +4h none 57.96821MS4 2s +2h 2s 2h 58.9841 MN4 5s +4h +p s +p 57.42383 MK4 2s +4h 2s 59.06624 Not included SN4 S4 SK4 Sixth Diurnals, total = 3, Species 6:M6 6s +6h none 86.952312MS6 4s +4h 2s 2h 87.96821 2MK6 4s +6h 2s 88.05035 Not included: MSN6 2SM6 MSK6 S6 2MN6 No others e.g. M8 Repeats:? Series6 ? 2MN2S2=ST36 Satellites:Coeff > 0.1 included 11/95(Doodson Numbers for p, N, p' to be added to astro argument of the major constituent, Cartwright coefficient applied to Amplitude of major constituent) 2Q1 0 1 0 0.1885 Sig1 0 1 0 0.1884 Q1 0 1 0 0.1884 RHO1 0 1 0 0.1882 O1 0 1 0 0.1885 NO1 0 1 0 0.2004 NO1 2 0 0 0.3596 S1 0 0 2 0.3534 K1 0 1 0 0.1356 J1 0 1 0 0.1356 OO1 2 0 0 0.1496 OO1 0 2 0 0.1342 OO1 0 1 0 0.6398 L2 2 1 0 +180 0.1102 L2 2 0 0 +180 0.2505 K2 0 1 0 0.2980 Satellites coeff > 0.035 or Latitude dependant NOT included 11/95 NO1 1 0 0 +90° 0.2227 Latitude dependant 2Q1 1 0 0 +270 0.0607 " " Rho1 2 0 0 +180 0.0576 J1 1 0 0 +90 0.0816 " " 2N2 0 1 0 +180 0.0374 2N2 1 0 0 +90 0.0678 " " Mu2 0 1 0 +180 0.0375 N2 0 1 0 +180 0.0373 Nu2 0 1 0 +180 0.0373 L2 0 1 0 +180 0.0366 M2 0 1 0 +180 0.0373 S2 0 +1 0 0.0811 SH; includes +2h! ++++++ Summary as used in .EXE DataBase ++++++ Note: Only upper case (capitals) used to facilitate programming Note: SGx numbers appear to end at 61 Format: [DATA] SGx# or 99etc , Constit Name , etc USED CONSTITUENTS: Total 37 DATA 1,O1,1,MK1,2,K1,2,MO1,2,SP1,3,M2,3,KO2,4,S2,4,KP2,5,M4,6,MS4,7,MF DATA 8,N2,8,KQ2,9,SA,10,MM,11,SIG1,11,SIGMA1,12,Q1,12,NK1,13,J1,13,MQ1 DATA 14,OO1,15,MU2,16,L2,16,2MN2,16,L2 2MN2,17,NO1,17,M1,18,MN4,20,SSA DATA 21,2Q1,22,RHO1 DATA 24,PI1,25,P1,25,SK1,26,S1,31,2N2,31,2MS2,32,NU2,34,T2,36,K2,37,M3 DATA 46,MO3,46,2MK3,47,SO3,48,MK3,51,MK4,55,M6,57,2MS6,58,2MK6,61,MSF UNUSED CONSTITUENTS: Total 108 4/12/95 99 or SGx used for known unused constits 98 used for total unknowns UK 97 used for total unknowns France DATA 19,M8,23,CHI1,27,PSI1,28,PHI1,29,THE1,30,OQ1,33,LDA2,35,R2,38,TAU1 DATA 39,SO1,40,EPS2,41,OP2,42,MKS2,43,MSN2,44,ETA2,45,2SM2,49,SK3,50,SN4 DATA 52,S4,53,SK4,54,2MN6,56,MSN6,59,2SM6,60,MSK6 end SGx numbers Total = 24 unused + 37 used = 61 DATA 98,MNUM,99,MSM Tot Long = 2 DATA 99,ALP1,99,UPS1,99,KQ1,99,BET1,99,MS1,99,MP1,98,TH1,97,KHI1,99,THETA1 Tot Diur = 9 DATA 99,OQ2,99,2MK2,99,SNK2,98,NA2,98,MA2,98,MB2,99,LAM2,99,2SK2,98,MSNU2 DATA 99,KJ2,99,SKM2,99,2MN2S2,99,2NS2,99,3M2S2,99,MNS2,99,MNUS2,99,M(SK)2 DATA 99,M(KS)2,99,LAMBDA2,99,NKM2 Tot Semis = 20 DATA 99,MQ3,98,2MS3,98,2MP3,99,MS3,98,2MQ3,99,SP3,99,S3,98,MP3,99,K3 Tot Thirds = 9 DATA 99,3MK4,99,3MS4,98,MSNK4,98,MNU4,99,2MSK4,98,MA4,99,2MKS4,99,3MN4 DATA 98,M2SK4,98,MR4,99,2MSN4,98,MT4,99,2NMS4,99,2MMUS4,99,2MNS4,99,2MNUS4 DATA 99,N4,97,3MN4ML4,97,NK4,99,2SMK4,99,2SNM4,99,2MKN4 Tot Quarters = 22 DATA 99,S6,99,4MK6,98,M2N6,99,4MS6,98,2MSNK6,98,2MNU6,98,3MSK6,98,MA6 DATA 98,3MKS6,98,4MN6,99,MNK6,98,2(MS)K6,98,2MT6,99,3MNS6,99,3MNUS6,99,2NM6 DATA 99,2SN6,99,3MSN6,99,3MKN6,99,3MNK6,97,4MN.2ML6,97,4MN.2ML Tot Sixths = 22 +++++++++++++++++++++++ NOTE 1 Sa :SH does not include p'; AD (in Admiralty Manual) does not include p' in Tables 7.1 or 15.1, but manages to reintroduce it in 7.17. MF appears to include it. Research problem. 20/12: include it.NOTE 2 M3 :Phase correction 180° or negative amplutude? Also AD implies only a quarter phase shift by addition of M1 and M2. 20/12: 180° in its own right. If interaction M1/M2, then should be found seperately.NOTE 3In general Cosine x° = Cosine x+360° = Cosine x360° (same holds for Sines). Therefore, is a phase shift of 90° equal to one of +270° ? Also is a phase shift of +/ 180° the equivalent of Amplitude * 1 ? I have a sneaking feeling at the back of my mind that phase shifts of 360° should not be admissable, but cannot prove it mathematically. 20/12: confirmed no diff.NOTE 4The only Latitude dependant satellite with coeff > 0.1 is third for NO1. Ask MF for details. 20/12: MF might be able to supply formulae.NOTE 5Astro: Cartwright appears to be in error with h (Lon Sun). Example: Table 3, p56, Geophys J R astr Soc 1971 gives Phi r line 3, col 2 as 60.11923°. Should be 61.38645° ???NOTE 6Stands: Is there a "standard definition" either in terms of time or height? Between well defined highs and lows, frequently diurnals, and as derivative (slope) approaches 0, with or without change of sign, should time be the only criterion?
