Harmonic Constituents for the free42
Overview
-"HCST" employs the gradient method to find the harmonic constituents
in a port, given 2000 water-levels.
-More than 2000 data may be used if you change the corresponding lines
in the listing.
Program Listing
-The water-levels may be calculated by:
WL(t) ~ SUM Aj [ cos ( fj
.t + Vj ) + aj cos ( fj .t + Vj
- 2 N' ) + bj cos ( fj .t + Vj -
N' ) + cj cos ( fj .t + Vj + N' ) +
dj cos ( fj .t + Vj + 2 N' ) ]
where a , b , c , d are constants to take the nodal corrections into account f = angular speed of the wave , V = phase
and N' = 234.98° + 0.05295376°
t where t = numbers of days since 2000/01/01 0h
UT
-You store in R100 thru R436 the approximate values of amplitudes A and phases V, and the gradient method improves these coefficients.
-I've used 2000 water-levels and 168 waves to get the harmonic constituents.
Data Registers: R00 = gradient-method parameter < 0 ( All the • Registers are to be initialized before executing "HCST" )
R01 = Sum errors^2 R02 = Sum errors R03 to R16: temp R17 thru R99: unused
-----------------------------------------------------------------------------------------------------------------------------------------------------
• R100 = Z0 • R101 = V1
• R103 = V2 .......................... •
R435 = V168
These registers are initialized by approximate values
• R102 = A1 • R104 = A2
.......................... • R436 = A168
and finally contain the improved constituents
------------------------------------------------------------------------------------------------------------------------------------------------------
• R501 = f1 • R506 = f2
.......................... • R1336 = f168
the 168 angular speeds of the 168 waves
• R502 = a1 • R507 = a2
.......................... • R1337 = a168
the 168 coefficients corresponding to -2.N'
• R503 = b1 • R508 = b2
.......................... • R1338 = b168
the 168 coefficients corresponding to -N'
• R504 = c1 • R509 = c2
.......................... • R1339 = c168
the 168 coefficients corresponding to +N'
• R505 = d1 • R510 = d2
.......................... • R1340 = d168
the 168 coefficients corresponding to +2.N'
------------------------------------------------------------------------------------------------------------------------------------------------------
• R1501 = t1 • R1502 = t2
.......................... • R3500 = t2000
the 2000 instants, expressed in days since 2000/01/01 0h UT
• R4001 = wl1 • R4002 = wl2 .......................
• R6000 = wl2000 the 2000 water-levels corresponding
to the instants above
R6501 = e1 R6502 = e2 ............................. R8500 = e2000 the 2000 errors
Flags: /
Subroutines: /
| 01 LBL "HCST"
02 -.0001 03 STO 00 04 234.98 05 STO 13 06 .05295376 07 STO 14 08 E99 09 STO 16 10 100.436 11 STO 11 12 LBL 01 13 8500 14 STO 03 15 6000 16 STO 04 17 3500 18 STO 05 19 CLX 20 STO 01 21 STO 02 22 LBL 02 23 RCL 11 24 STO 06 |
25 501
26 STO 07 27 RCL IND 06 28 STO 08 29 ISG 06 30 LBL 03 31 RCL IND 05 32 RCL IND 07 33 * 34 RCL IND 06 35 + 36 STO 09 37 COS 38 LASTX 39 RCL 14 40 RCL IND 05 41 * 42 RCL 13 43 + 44 STO 10 45 ST+ X 46 - 47 COS 48 ISG 07 |
49 CLX
50 RCL IND 07 51 * 52 + 53 RCL 09 54 RCL 10 55 - 56 COS 57 ISG 07 58 CLX 59 RCL IND 07 60 * 61 + 62 RCL 09 63 RCL 10 64 + 65 COS 66 ISG 07 67 CLX 68 RCL IND 07 69 * 70 + 71 RCL 09 72 RCL 10 |
73 ST+ X
74 + 75 COS 76 ISG 07 77 CLX 78 RCL IND 07 79 * 80 + 81 ISG 06 82 RCL IND 06 83 * 84 ST+ 08 85 ISG 07 86 CLX 87 ISG 06 88 GTO 03 89 RCL 08 90 RCL IND 04 91 - 92 STO IND 03 93 ST+ 02 94 X^2 95 ST+ 01 96 DSE 03 |
97 DSE 04
98 DSE 05 99 4000 100 RCL 04 101 X>Y? 102 GTO 02 103 RCL 00 104 RCL 02 105 * 106 ST+ IND 11 107 RCL 11 108 STO 06 109 501 110 STO 07 111 ISG 06 112 LBL 04 113 CLX 114 STO 12 115 STO 15 116 3500 117 STO 05 118 8500 119 STO 03 120 LBL 05 |
121 RCL IND 05
122 RCL IND 07 123 * 124 RCL IND 03 125 P-R 126 ST+ 12 127 X<>Y 128 ST+ 15 129 DSE 03 130 DSE 05 131 1500 132 RCL 05 133 X>Y? 134 GTO 05 135 RCL 00 136 ST* 12 137 ST* 15 138 RCL IND 06 139 ISG 06 140 RCL IND 06 141 P-R 142 X<>Y 143 RCL 15 144 - |
145 X<>Y
146 RCL 12 147 + 148 R-P 149 STO IND 06 150 X<>Y 151 DSE 06 152 STO X 153 STO IND 06 154 CLX 155 SIGN 156 ST+ 06 157 5 158 ST+ 07 159 ISG 06 160 GTO 04 161 RCL 01 162 ENTER 163 X<> 16 164 X<>Y 165 VIEW 01 166 X<Y? 167 GTO 01 168 END |
( 285 bytes / SIZE 8501 )
| STACK | INPUTS | OUTPUTS |
| X | / | SUM errors^2 |
Example: Find the harmonic
constituents in Brest ( France )
-The French "SHOM" ( Service Hydrographique et Océanographique
de la Marine ) gives the water-levels in many ports.
-Unfortunately - for an unknown reason - they refuse to publish the
constants and even the 143 waves that they employ.
-So, I've used 168 waves - listed below - to estimate the harmonic constituents from 2000 water-levels in 2018-2019:
-With d = the number of days since 01/01/2000 at 0h UT
h = (279.97°) + 0.98564736°
d
( all angles are expressed in degrees )
s = (211.73°) + 13.17639648° d
p = (83.30°) + 0.11140352°
d
p1 = (282.94°) + 0.00004708°
d
With 168 WAVES ( 7+21+39+12+23+10+22+5+16+9+4 )
Wave Speed (deg/day) a/b/c/d Expression
SA 0.98564736
0/0/0/0
h
SSA 1.97129472
0/0/0/0
2h
MSM 11.31650528 0/-0.072/-0.065/0
s-2h+p
Mm 13.06499296 0/-0.066/-0.065/0
s-p
MSf 24.38149824 0/-0.072/-0.065/0
2s-2h
Mf 26.35279296
0/0/0.415/0.039 2s
MTM 39.41778592 0/0/0.414/0.039
3s-p
-------------------------------------------------------
2Q1 308.50286848 0/0.189/0/0
360d-4s+h+2p
SIG1 310.25135616 -0.006/0.189/0/0
360d-4s+3h
Q1 321.56786144 -0.006/0.189/0/0
360d-3s+h+p
RHO1 323.31634912 -0.005/0.189/0/0
360d-3s+3h-p
O1 334.63285440 -0.006/0.189/0/0
360d-2s+h
MP1 336.60414912 0/-0.037/-0.011/0
360d-2s+3h
A19 345.94935968 -0.013/-0.229/0/0
360d-s-h+p
M1 347.92065440 0/-0.029/0.201/-0.006
360d-s+h+p
M1B 347.69784736 -0.016/0.185/0/0
360d-s+h-p
CHI1 349.66914208 0/-0.029/0.219/-0.005 360d-s+3h-p
PI1 358.02875236 0/0/0/0
360d-2s+p1
P1 359.01435264 0/0/0/0
360d-h
S1 360
0/0/0/0
360d
K1 360.98564736 0/-0.020/0.136/-1:343
360d+h
PSI1 361.97124764 0/0.004/0.019/0
360d+2h-p1
PHI1 362.95694208 0/0/-0.039/-0.019
360d+3h
THET1 372.30215264 0/-0.032/0.198/-0.004 360d+s-h+p
J1 374.05064032 0/-0.029/0.199/0
360d+s+h-p
SO1 385.36714560 0/0/0.189/0
360d+2s-h
OO1 387.33844032 0/0/0.640/0.134/0.009 360d+2s+h
UPS1 400.40343328 0/0/0.640/0.134
360d+3s+h-p
--------------------------------------------------------
2MN2S2 633.79051232 0/-0.112/0.004/0
720d-7s+6h+p
3MKS2 644.88421056 0.032/0.186/-0.011/0 720d-6s+4h
2NS2 645.10701760 0/-0.075/0/0
720d-6s+4h+2p
3M2S2 646.85550528 0/-0.112/0.004/0
720d-6s+6h
OQ2 656.20071584 -0.012/0.378/0/0
720d-5s+2h+p
??? 656.42352208
0/0/0/0
720d-5s+2h+3p
MNS2 658.17201056 0/-0.075/0/0
720d-5s+4h+p
MNUS2 659.92049824 0/-0.075/0.002/0
720d-5s+6h-p
2MS2K2 667.29441408 0.032/0.596/-0.026/0 720d-4s
2MK2 669.26570880 0.032/0.223/-0.013/0
720d-4s+2h
2N2 669.48851584
0/-0.037/0/0
720d-4s+2h+2p
MU2 671.23700352 0/-0.038/0/0
720d-4s+4h
??? 672.22260380
0/-0.037/0/0
720d-4s+5h-p1
NA2 681.56790852 0/-0.032/0/0
720d-3s+h+p+p1
N2 682.55350880
0/-0.037/0/0
720d-3s+2h+p
NU2 684.30199648 0/-0.038/0/0
720d-3s+4h-p
MKL2S2 686.27329128 0/-0.088/0.302/0.032 720d-3s+6h-p
OP2 693.64720704 0/0.178/0/0
720d-2s
GAM2 693.87001408 0/-0.035/0/0
720d-2s+2p
ALPH2 694.63290148 0/0/0/0
720d-2s+h+p1
M2 695.61850176
0/-0.037/0/0
720d-2s+2h
BET2 696.60410204 0/-0.021/0/0
720d-2s+3h-p1
DEL2 697.58979648 0/-0.037/0/0
720d-2s+4h
M2KS2 699.56109120 0/-0.063/0.600/0.064
720d-2s+6h
2SNMK2 704.96371232 0.032/0.227/-0.050/0 720d-s-2h+p
LAM2 706.93500704 0/-0.045/0/0
720d-s+p
L2 708.68349472
0/-0.037/0/0
720d-s+2h-p
L2B 708.90630172 0/-0.015/0.441/0.062
720d-s+2h+p
T2 719.01439972
0/0/0/0
720d-h+p1
S2
720 0/0.002/0/0
720d
R2 720.98560028
0/0/0/0
720d+h-p1
K2 721.97129472
0/-0.013/0.298/0.032 720d+2h
MSNU2 731.31650528 0/-0.035/-0.038/0
720d+s-2h+p
MSN2 733.06499296 0/-0.035/-0.037/0
720d+s-p
XI2 733.28780000
0/-0.021/0.436/0.047 720d+s+p
KJ2 735.03628768
0/-0.049/0.335/-0.003 720d+s+2h-p
2KMSN2 737.00758240 0/-0.063/0.561/0.064
720d+s+4h-p
2SM2 744.38149824 0/-0.004/-0.037/0
720d+2s-2h
SKM2 746.35279296 0/-0.011/0.261/0.032
720d+2s
------------------------------------------------------------
MQ3 1017.18636320 -0.006/0.152/0/0
1080d-5s+3h+p
2MK3 1030.25135616 0/0.061/-0.020/0
1080d-4s+3h
2MP3 1032.22265088 0/-0.075/0/0
1080d-4s+5h
M3 1043.42775264 0/-0.056/0/0
1080d-3s+3h
SO3 1054.63285440 -0.006/0.191/0/0
1080d-2s+h
MK3 1056.60414912 0/-0.057/0.136/0
1080d-2s+3h
A87 1069.78054560 0/-0.038/0.439/0.048
1080d-s+3h
SP3 1079.01435264 0/-0.009/0/0
1080d-h
S3 1080
0/0/0/0
1080d
SK3 1080.98564736 0/-0.018/0.136/-1:343
1080d+h
K3 1082.95694208 0/-0.060/0.408/0.010
1080d+3h
2SO3 1105.36714560 0/0.004/0.189/-0.006
1080d+2s-h
------------------------------------------------------------
2NMS4 1340.72551936 0/-0.112/0.002/0
1440d-8s+6h+2p
2MMUS4 1342.47400704 0/-0.112/0.002/0
1440d-8s+8h
2MNS4 1353.79051232 0/-0.112/0.002/0
1440d-7s+6h+p
2MNUS4 1355.53900000 0/-0.112/0.002/0
1440d-7s+8h-p
3MK4 1364.88421056 0.032/0.187/-0.013/0
1440d-6s+4h
N4 1365.10701760 0/-0.075/0/0
1440d-6s+4h+2p
3MS4 1366.85550528 0/-0.112/0.002/0
1440d-6s+6h
MN4 1378.17201056 0/-0.075/0/0
1440d-5s+4h+p
MNU4 1379.92049824 0/-0.075/0/0
1440d-5s+6h-p
2MSK4 1389.26570880 0.032/0.223/-0.013/0
1440d-4s+2h
M4 1391.23700352 0/-0.075/0/0
1440d-4s+4h
2MKS4 1393.20829824 0/-0.088/0.298/0.002
1440d-4s+6h
SN4 1402.55350880 0/-0.035/0/0
1440d-3s+2h+p
3MN4 1404.30199648 0/-0.112/-0.037/0
1440d-3s+4h-p
NK4 1404.52480352 0/-0.050/0.298/0
1440d-3s+4h+p
MT4 1414.63290148 0/-0.037/0/0
1440d-2s+h+p1
MS4 1415.61850176 0/-0.035/0/0
1440d-2s+2h
MK4 1417.58979648 0/-0.050/0.298/0
1440d-2s+4h
2SNM4 1426.93500704 0/-0.035/-0.037/0
1440d-s+p
2MSN4 1428.68349472 0/-0.075/-0.037/0
1440d-s+2h-p
2MKN4 1430.65478944 0/-0.088/0.261/0.032
1440d-s+4h-p
S4 1440
0/0/0/0
1440d
SK4 1441.97129472 0/-0.011/0.298/0.032
1440d+2h
--------------------------------------------------------------
MNO5 1712.80486496 0/0.115/0/0
1800d-7s+5h+p
2MO5 1725.86985792 0/0.115/0/0
1800d-6s+5h
3MP5 1727.84115264 0/-0.112/-0.011/0
1800d-6s+7h
M5 1739.04625440 0/-0.094/0/0
1800d-5s+5h
MNK5 1739.15765792 0/-0.095/0.134/-0.003
1800d-5s+5h+p
2MP5 1750.25135616 0/-0.086/0/0
1800d-4s+3h
2MK5 1752.22265088 0/-0.095/0.136/-0.003
1800d-4s+5h
MSK5 1776.60414912 0/-0.055/0.136/-0.003
1800d-2s+3h
3KM5 1778.57544384 0/-0.097/0.408/-0.009
1800d-2s+5h
2SK5 1800.98564736 0/-0.016/0.136/-0.003
1800d+h
--------------------------------------------------------------
3MNK6 2047.43771936 0.032/0.186/-0.013/0
2160d-9s+6h+p
3MNS6 2049.40901408 0/-0.149/0.002/0
2160d-9s+8h+p
3MNUS6 2051.15750176 0/-0.150/0.002/0
2160d-9s+10h-p
4MK6 2060.50271232 0.032/0.149/-0.013/0
2160d-8s+6h
2NM6 2060.72551936 0/-0.112/0/0
2160d-8s+6h+2p
4MS6 2062.47400704 0/-0.147/0/0
2160d-8s+8h
2MN6 2073.79051232 0/-0.112/0/0
2160d-7s+6h+p
2MNU6 2075.53900000 0/-0.113/0/0
2160d-7s+8h-p
3MSK6 2084.88421056 0.032/0.188/-0.013/0
2160d-6s+4h
M6 2086.85550528 0/-0.112/0/0
2160d-6s+6h
3MKS6 2088.82680000 0/-0.125/0.300/0.032
2160d-6s+8h
MSN6 2098.17201056 0/-0.073/0/0
2160d-5s+4h+p
4MN6 2099.92049824 0/-0.149/-0.037/0
2160d-5s+6h-p
MNK6 2100.14330528 0/-0.088/0.298/0.032
2160d-5s+6h+p
2MT6 2110.25140324 0/-0.075/0/0
2160d-4s+3h+p1
2MS6 2111.23700352 0/-0.073/0/0
2160d-4s+4h
2MK6 2113.20829824 0/-0.088/0.298/0.032
2160d-4s+6h
2SN6 2122.55350880 0/-0.033/0/0
2160d-3s+2h+p
3MSN6 2124.30199648 0/-0.110/-0.037/0
2160d-3s+4h-p
3MKN6 2126.27329120 0/-0.125/0.261/0.032
2160d-3s+6h-p
2SM6 2135.61850176 0/-0.033/0/0
2160d-2s+2h
MSK6 2137.58979648 0/-0.048/0.298/0.032
2160d-2s+4h
--------------------------------------------------------------
2MNO7 2408.42336672 -0.006/0.077/0/0
2520d-9s+7h+p
2NMK7 2421.71116672 0/-0.132/0.136/-0.003
2520d-8s+7h+2p
M7 2434.77615968 0/-0.129/0/0
2520d-7s+7h+p
2MSO7 2445.86985792 -0.006/0.116/0/0
2520d-6s+5h
MSKO7 2472.22265088 -0.006/0.141/0.298/0.032 2520d-4s+5h
--------------------------------------------------------------
2MNS8 2745.02751584 0/-0.110/0/0
2880d-11s+10h+p
2(MN)8 2756.34402112 0/-0.149/0/0
2880d-10s+8h+2p
3MN8 2769.40901408 0/-0.149/0/0
2880d-9s+8h+p
3MNU8 2771.15750176 0/-0.149/0/0
2880d-9s+10h-p
4MSK8 2780.50271232 0.032/0.151/-0.013/0
2880d-8s+6h
M8 2782.47400704 0/-0.149/0/0
2880d-8s+8h
2MSN8 2793.79051232 0/-0.109/0/0
2880d-7s+6h+p
3ML8 2795.53900000 0/-0.149/0/0
2880d-7s+8h-p
2MNK8 2795.76180704 0/-0.125/0.298/0.032
2880d-7s+8h+p
3MS8 2806.85550528 0/-0.110/0/0
2880d-6s+6h
3MK8 2808.82680000 0/-0.125/0.298/0.032
2880d-6s+8h
2SMN8 2818.17201056 0/-0.071/0/0
2880d-5s+4h+p
4MSN8 2819.92049824 0/-0.147/-0.037/0
2880d-5s+6h-p
MSNK8 2820.14330528 0/-0.086/0.298/0.032
2880d-5s+6h+p
2(MS)8 2831.23700352 0/-0.071/0/0
2880d-4s+4h
2MSK8 2833.20829824 0/-0.086/0.298/0.032
2880d-4s+6h
---------------------------------------------------------------
5MNS10 3440.64601760 0/-0.224/0.002/0
3600d-13s+12h+p
3M2N10 3451.96251288 0/-0.187/0/0
3600d-12s+10h+2p
4MN10 3465.02751584 0/-0.187/0/0
3600d-11s+10h+p
M10 3478.09250880 0/-0.187/0/0
3600d-10s+10h
3MSN10 3489.40901408 0/-0.147/0/0
3600d-9s+8h+p
4MS10 3502.47400704 0/-0.147/0/0
3600d-8s+8h
4MK10 3504.44530176 0/-0.162/0.298/0.032
3600d-8s+10h
5MSN10 3515.53900000 0/-0.186/-0.037/0
3600d-7s+8h-P
3M2S10 3526.85550528 0/-0.098/0/0
3600d-6s+6h
---------------------------------------------------------------
4MNS12 4185.02751584 0/-0.183/0/0
4320d-11s+10h+p
5MS12 4198.09250880 0/-0.183/0/0
4320d-10s+10h
4MSL12 4211.15750176 0/-0.183/0/0
4320d-9s+10h-p
4M2S12 4222.47400704 0/-0.146/0/0
4320d-8s+8h
-Store the angular speeds of the waves and the coefficients a b c d
in registers R501 to R1340
-Then, the times ( in days since 2000/01/01 0h UT ) in R1501 thru R3500
and the corresponding water-levels in R4001 thru R6000
-If you have approximate values for the phases and amplitudes, store
them in R101 thru R436
( otherwise, store an approximate value of the mean-sea-level
in R100 & zero in R101 to R436 )
-Then, XEQ "HCST" >>>> free42 will display the successive R01
= sum of errors^2
-The program will stop when R01 starts to increase ( it may last 1
or 2 hours )
-If the routine doesn't stop, stop it when R01 becomes almost constant
and enough small.
-You can also set flag F21 and "HCST" will stop at each iteration.
-It will work with free42decimal but free42binary is faster !
-With 168 waves, each iteration lasts about 5 seconds with free42binary.
-The gradient method employs the 1st derivatives. I've neglected the
nodal corrections a b c d to compute them.
-For the "gradient parameter" in R00, I've used -0.0001 but
another value is perhaps better:
-Change line 02 or store another value in R00 before running the routine
again.
>>>> The results I've obtained are used in the routine "BREST" listed
in "Brest Water-Level for the HP41"
References:
[1] Bernard Simon - " La Marée Océanique Côtière"
[2] Harmonic
Constituents with Nodal Corrections - November 2006 - Edition n°1
[3] Hartmann and Wenzel - ( 1995 ) - Tidal Potential Catalogue ( HW95
)
[4] Cartwright, Tayler and Edden - Tidal Potential Catalogue
( 1971 & 1973 ).