1   /* Copyright 2002-2024 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.models.earth.displacement;
18  
19  import org.hipparchus.exception.LocalizedCoreFormats;
20  import org.hipparchus.util.FastMath;
21  import org.orekit.bodies.GeodeticPoint;
22  import org.orekit.data.DataSource;
23  import org.orekit.errors.OrekitException;
24  import org.orekit.errors.OrekitMessages;
25  
26  import java.io.BufferedReader;
27  import java.io.IOException;
28  import java.util.ArrayList;
29  import java.util.List;
30  import java.util.regex.Matcher;
31  import java.util.regex.Pattern;
32  
33  /**
34   * Parser for ocean loading coefficients, using Onsala Space Observatory files in BLQ format.
35   * <p>
36   * Files in BLQ format can be generated using the form at the
37   * <a href="http://holt.oso.chalmers.se/loading/">Bos-Scherneck web site</a>,
38   * selecting BLQ as the output format.
39   * </p>
40   * <p>
41   * The sites names are extracted from the file content, not the file name, because the
42   * file can contain more than one station. As we expect existing files may have been
43   * stripped from headers and footers, we do not attempt to parse them. We only parse
44   * the series of 7 lines blocks starting with the lines with the station names and their
45   * coordinates and the 6 data lines that follows. Several such blocks may appear in the
46   * file. Copy-pasting the entire mail received from OSO after completing the web site
47   * form works, as intermediate lines between the 7 lines blocks are simply ignored.
48   * </p>
49   * @see OceanLoadingCoefficients
50   * @see OceanLoading
51   * @since 9.1
52   * @author Luc Maisonobe
53   */
54  public class OceanLoadingCoefficientsBlqParser {
55  
56      /** Pattern for fields with real type. */
57      private static final String  REAL_TYPE_PATTERN = "[-+]?(?:\\p{Digit}+(?:\\.\\p{Digit}*)?|\\.\\p{Digit}+)(?:[eE][-+]?\\p{Digit}+)?";
58  
59      /** Pattern for extracted real fields. */
60      private static final String  REAL_FIELD_PATTERN = "\\p{Space}*(" + REAL_TYPE_PATTERN + ")";
61  
62      /** Pattern for end of line. */
63      private static final String  END_OF_LINE_PATTERN = "\\p{Space}*$";
64  
65      /** Pattern for site name and coordinates lines. */
66      private static final String  SITE_LINE_PATTERN = "^\\$\\$ *([^,]*),\\p{Space}*(?:RADI TANG)?\\p{Space}*lon/lat:" +
67                                                       REAL_FIELD_PATTERN +
68                                                       REAL_FIELD_PATTERN +
69                                                       REAL_FIELD_PATTERN +
70                                                       END_OF_LINE_PATTERN;
71  
72      /** Pattern for coefficients lines. */
73      private static final String  DATA_LINE_PATTERN = "^" +
74                                                       REAL_FIELD_PATTERN + // M₂ tide
75                                                       REAL_FIELD_PATTERN + // S₂ tide
76                                                       REAL_FIELD_PATTERN + // N₂ tide
77                                                       REAL_FIELD_PATTERN + // K₂ tide
78                                                       REAL_FIELD_PATTERN + // K₁ tide
79                                                       REAL_FIELD_PATTERN + // O₁ tide
80                                                       REAL_FIELD_PATTERN + // P₁ tide
81                                                       REAL_FIELD_PATTERN + // Q₁ tide
82                                                       REAL_FIELD_PATTERN + // Mf tide
83                                                       REAL_FIELD_PATTERN + // Mm tide
84                                                       REAL_FIELD_PATTERN + // Ssa tide
85                                                       END_OF_LINE_PATTERN;
86  
87      /** Pattern for site name and coordinates lines. */
88      private static final Pattern SITE_PATTERN = Pattern.compile(SITE_LINE_PATTERN);
89  
90      /** Pattern for coefficients lines. */
91      private static final Pattern DATA_PATTERN = Pattern.compile(DATA_LINE_PATTERN);
92  
93      /** Main tides. */
94      private static final Tide[][] TIDES = {
95          {
96              Tide.SSA, Tide.MM, Tide.MF
97          }, {
98              Tide.Q1,  Tide.O1, Tide.P1, Tide.K1
99          }, {
100             Tide.N2,  Tide.M2, Tide.S2, Tide.K2
101         }
102     };
103 
104     /** Species index for each column. */
105     private static final int[] I = {
106         2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0
107     };
108 
109     /** Tides index for each column. */
110     private static final int[] J = {
111         1, 2, 0, 3, 3, 1, 2, 0, 2, 1, 0
112     };
113 
114     /** Parse a BLQ file.
115      * <p>
116      * Files in BLQ format can be generated using the form at the
117      * <a href="https://holt.oso.chalmers.se/loading/">Bos-Scherneck web site</a>,
118      * selecting BLQ as the output format.
119      * </p>
120      * <p>
121      * when completing the web site form, the email received as the following form:
122      * </p>
123      * <pre>{@literal
124      * $$ Ocean loading displacement
125      * $$
126      * $$ Calculated on holt using olfg/olmpp of H.-G. Scherneck
127      * $$
128      * $$ COLUMN ORDER:  M2  S2  N2  K2  K1  O1  P1  Q1  MF  MM SSA
129      * $$
130      * $$ ROW ORDER:
131      * $$ AMPLITUDES (m)
132      * $$   RADIAL
133      * $$   TANGENTL    EW
134      * $$   TANGENTL    NS
135      * $$ PHASES (degrees)
136      * $$   RADIAL
137      * $$   TANGENTL    EW
138      * $$   TANGENTL    NS
139      * $$
140      * $$ Displacement is defined positive in upwards, South and West direction.
141      * $$ The phase lag is relative to Greenwich and lags positive. The
142      * $$ Gutenberg-Bullen Greens function is used. In the ocean tide model the
143      * $$ deficit of tidal water mass has been corrected by subtracting a uniform
144      * $$ layer of water with a certain phase lag globally.
145      * $$
146      * $$ Complete <model name> : No interpolation of ocean model was necessary
147      * $$ <model name>_PP       : Ocean model has been interpolated near the station
148      * $$                         (PP = Post-Processing)
149      * $$
150      * $$ Ocean tide model: CSR4.0, long-period tides from FES99
151      * $$
152      * $$ END HEADER
153      * $$
154      *   Goldstone
155      * $$ Complete CSR4.0_f
156      * $$ Computed by OLFG, H.-G. Scherneck, Onsala Space Observatory 2017-Sep-28
157      * $$ Goldstone,                 RADI TANG  lon/lat:  243.1105   35.4259    0.000
158      *   .00130 .00155 .00064 .00052 .01031 .00661 .00339 .00119 .00005 .00002 .00003
159      *   .00136 .00020 .00024 .00004 .00322 .00202 .00106 .00036 .00007 .00003 .00001
160      *   .00372 .00165 .00082 .00045 .00175 .00113 .00057 .00022 .00004 .00002 .00003
161      *     -2.7 -106.3  -62.6 -106.8   41.6   27.3   40.4   24.0 -119.1 -123.2 -169.7
162      *   -145.3  -88.4  178.5  -66.3 -130.5 -145.3 -131.7 -148.7  124.3  139.6   23.3
163      *     90.7  111.1   74.1  111.3  176.9  165.3  175.8  164.0   48.9   25.3    4.5
164      * $$
165      *   ONSALA
166      * $$ CSR4.0_f_PP ID: 2017-09-28 15:01:14
167      * $$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2017
168      * $$ Onsala,                    RADI TANG  lon/lat:   11.9264   57.3958    0.000
169      *   .00344 .00121 .00078 .00031 .00189 .00116 .00064 .00004 .00090 .00048 .00041
170      *   .00143 .00035 .00035 .00008 .00053 .00051 .00018 .00009 .00013 .00006 .00007
171      *   .00086 .00023 .00023 .00006 .00029 .00025 .00010 .00008 .00003 .00001 .00000
172      *    -64.6  -50.3  -95.0  -53.1  -58.8 -152.4  -65.5 -133.8    9.8    5.8    2.1
173      *     85.4  115.2   56.7  114.7   99.5   15.9   94.2  -10.0 -166.3 -169.8 -177.7
174      *    110.7  147.1   93.9  148.6   49.4  -56.5   34.8 -169.9  -35.3   -3.7   10.1
175      * $$ END TABLE
176      * Errors:
177      * Warnings:
178      * }</pre>
179      * <p>
180      * We only parse blocks 7 lines blocks starting with the lines with the station names
181      * and their coordinates and the 6 data lines that follows. Several such blocks may
182      * appear in the file.
183      * </p>
184      * @param source source for BLQ data
185      * @return parsed coefficients
186      */
187     public List<OceanLoadingCoefficients> parse(final DataSource source) {
188 
189         final List<OceanLoadingCoefficients> coefficients = new ArrayList<>();
190 
191         // temporary holders for parsed data
192         String siteName = null;
193         GeodeticPoint siteLocation = null;
194         final double[][][] data = new double[6][3][];
195         for (int i = 0; i < data.length; ++i) {
196             data[i][0] = new double[3];
197             data[i][1] = new double[4];
198             data[i][2] = new double[4];
199         }
200 
201         try (BufferedReader reader = new BufferedReader(
202             source.getOpener().openReaderOnce())) {
203 
204             int lineNumber = 0;
205             int dataLine = -1;
206             for (String line = reader.readLine(); line != null; line = reader.readLine()) {
207                 ++lineNumber;
208 
209                 if (dataLine < 0) {
210                     // we are looking for a site line
211                     final Matcher matcher = SITE_PATTERN.matcher(line);
212                     if (matcher.matches()) {
213                         // the current line is a site description line
214                         siteName = matcher.group(1);
215                         siteLocation =
216                             new GeodeticPoint(FastMath.toRadians(Double.parseDouble(matcher.group(3))),
217                                               FastMath.toRadians(Double.parseDouble(matcher.group(2))),
218                                               Double.parseDouble(matcher.group(4)));
219                         // next line must be line 0 of the data
220                         dataLine = 0;
221                     }
222                 } else {
223                     // we are looking for a data line
224                     final Matcher matcher = DATA_PATTERN.matcher(line);
225                     if (!matcher.matches()) {
226                         throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
227                                                   lineNumber, source.getName(), line);
228                     }
229                     for (int k = 0; k < I.length; ++k) {
230                         if (dataLine < 3) {
231                             // amplitude data
232                             data[dataLine][I[k]][J[k]] = Double.parseDouble(matcher.group(k + 1));
233                         } else {
234                             // phase data (reversed to be negative for lags)
235                             data[dataLine][I[k]][J[k]] = -FastMath.toRadians(Double.parseDouble(matcher.group(k + 1)));
236                         }
237                     }
238                     if (dataLine < data.length - 1) {
239                         // we need more data
240                         ++dataLine;
241                     } else {
242                         // it was the last data line
243                         coefficients.add(new OceanLoadingCoefficients(siteName, siteLocation,
244                                                                       TIDES, data[0],
245                                                                       data[3], data[1],
246                                                                       data[4], data[2],
247                                                                       data[5]));
248                         dataLine = -1;
249                     }
250                 }
251             }
252 
253             if (dataLine >= 0) {
254                 // we were looking for a line that did not appear
255                 throw new OrekitException(OrekitMessages.UNEXPECTED_END_OF_FILE_AFTER_LINE,
256                                           source.getName(), lineNumber);
257             }
258 
259         } catch (IOException ioe) {
260             throw new OrekitException(ioe, LocalizedCoreFormats.SIMPLE_MESSAGE,
261                                       ioe.getLocalizedMessage());
262         }
263 
264         return coefficients;
265 
266     }
267 
268 }