001/* 002 * This file is part of McIDAS-V 003 * 004 * Copyright 2007-2024 005 * Space Science and Engineering Center (SSEC) 006 * University of Wisconsin - Madison 007 * 1225 W. Dayton Street, Madison, WI 53706, USA 008 * https://www.ssec.wisc.edu/mcidas/ 009 * 010 * All Rights Reserved 011 * 012 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and 013 * some McIDAS-V source code is based on IDV and VisAD source code. 014 * 015 * McIDAS-V is free software; you can redistribute it and/or modify 016 * it under the terms of the GNU Lesser Public License as published by 017 * the Free Software Foundation; either version 3 of the License, or 018 * (at your option) any later version. 019 * 020 * McIDAS-V is distributed in the hope that it will be useful, 021 * but WITHOUT ANY WARRANTY; without even the implied warranty of 022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 023 * GNU Lesser Public License for more details. 024 * 025 * You should have received a copy of the GNU Lesser Public License 026 * along with this program. If not, see https://www.gnu.org/licenses/. 027 */ 028 029package edu.wisc.ssec.mcidasv.adt; 030 031import static java.lang.Math.asin; 032import static java.lang.Math.sin; 033 034import java.util.Scanner; 035 036public class Functions { 037 038 static String[] Months = { "JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", 039 "OCT", "NOV", "DEC" }; 040 041 static double[] PW_TnoValues = { -9999., -8888., 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 042 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 043 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 044 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 045 7.0, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 046 8.7, 8.8, 8.9, 9.0 }; 047 048 static double[][] PW_PressureValues = { 049 /* Atlantic pressure relationship values */ 050 { -9999.0, -8888.0, 1014.0, 1013.6, 1013.2, 1012.8, 1012.4, 1012.0, 1011.4, 1010.8, 051 1010.2, 1009.6, 1009.0, 1008.2, 1007.4, 1006.6, 1005.8, 1005.0, 1004.0, 1003.0, 052 1002.0, 1001.0, 1000.0, 998.8, 997.6, 996.4, 995.2, 994.0, 992.6, 991.2, 989.8, 053 988.4, 987.0, 985.4, 983.8, 982.2, 980.6, 979.0, 977.2, 975.4, 973.6, 971.8, 054 970.0, 968.0, 966.0, 964.0, 962.0, 960.0, 957.6, 955.2, 952.8, 950.4, 948.0, 055 945.4, 942.8, 940.2, 937.6, 935.0, 932.2, 929.4, 926.6, 923.8, 921.0, 918.0, 056 915.0, 912.0, 909.0, 906.0, 902.8, 899.6, 896.4, 893.2, 890.0, 886.6, 883.2, 057 879.8, 876.4, 873.0, 869.4, 865.8, 862.2, 858.6, 855.0 }, 058 /* Pacific pressure relationship values */ 059 { -9999.0, -8888.0, 1005.0, 1004.6, 1004.2, 1003.8, 1003.4, 1003.0, 1002.4, 1001.8, 060 1001.2, 1000.6, 1000.0, 999.4, 998.8, 998.2, 997.6, 997.0, 995.8, 994.6, 993.4, 061 992.2, 991.0, 989.6, 988.2, 986.8, 985.4, 984.0, 982.4, 980.8, 979.2, 977.6, 062 976.0, 974.0, 972.0, 970.0, 968.0, 966.0, 963.6, 961.2, 958.8, 956.4, 954.0, 063 951.4, 948.8, 946.2, 943.6, 941.0, 938.2, 935.4, 932.6, 929.8, 927.0, 924.4, 064 921.8, 919.2, 916.6, 914.0, 910.8, 907.6, 904.4, 901.2, 898.0, 894.2, 890.4, 065 886.6, 882.8, 879.0, 874.8, 870.6, 866.4, 862.2, 858.0, 853.4, 848.8, 844.2, 066 839.6, 835.0, 830.0, 825.0, 820.0, 815.0, 810.0 } }; 067 068 /** Atlantic/Pacific pressure relationship values */ 069 static double[] PW_WindValues = { -9999.0, -8888.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 26.0, 070 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0, 35.0, 37.0, 39.0, 41.0, 43.0, 45.0, 071 47.0, 49.0, 51.0, 53.0, 55.0, 57.0, 59.0, 61.0, 63.0, 65.0, 67.4, 69.8, 72.2, 74.6, 072 77.0, 79.6, 82.2, 84.8, 87.4, 90.0, 92.4, 94.8, 97.2, 99.6, 102.0, 104.6, 107.2, 109.8, 073 112.4, 115.0, 117.4, 119.8, 122.2, 124.6, 127.0, 129.6, 132.2, 134.8, 137.4, 140.0, 074 143.0, 146.0, 149.0, 152.0, 155.0, 158.0, 161.0, 164.0, 167.0, 170.0, 173.0, 176.0, 075 179.0, 182.0, 185.0, 188.0, 191.0, 194.0, 197.0, 200.0 }; 076 077 private static double RADIANSCONSTANT = 0.017453292; 078 079 private static double ANGLE90DEGRADIANS = 1.570797; 080 081 private static double EARTHRADIUSKM = 6371.0; 082 083 private static int[] JulianDateMonthArray = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 084 334, 365 }; 085 086 private static double deg_to_rad(double deg) { 087 return ((deg * Math.PI) / 180.0); 088 } 089 090 private static double rad_to_deg(double rad) { 091 return ((rad * 180.0) / Math.PI); 092 } 093 094 private static double a_sign(double x, double y) { 095 return ((x) / (y)) * Math.abs(y); 096 } 097 098 /** 099 * Compute time in ADT xxxxx.yyy units, where xxxxx is the day and yyy is 100 * the percentage of the day. 101 * 102 * This routine will also correct for Y2K problems. 103 * 104 * @param InputJulianDate 105 * Julian date. 106 * @param InputHMSTime 107 * Time in HHMMSS format. 108 * 109 * @return Date/Time value in xxxxx.yyy units. 110 */ 111 public static double calctime(int InputJulianDate, int InputHMSTime) { 112 double TimeReturnValue; 113 114 /* System.out.printf("calctime : input date=%d\n", InputJulianDate); */ 115 if (((InputJulianDate % 1000) == 0) || (InputHMSTime < 0)) { 116 TimeReturnValue = 0.0; 117 } else { 118 int YearValue = InputJulianDate / 1000; 119 /** obtain year */ 120 /* 121 * check for millenium designation in the year. if it is not there, 122 * add it onto the beginning 123 */ 124 /* System.out.printf("calctime : year=%d\n", YearValue); */ 125 if (YearValue < 1900) { 126 if (YearValue > 70) { 127 InputJulianDate = 1900000 + InputJulianDate; 128 } else { 129 InputJulianDate = 2000000 + InputJulianDate; 130 } 131 } 132 /* 133 * System.out.printf("calctime : InputJulianDate=%d\n",InputJulianDate 134 * ); 135 */ 136 double SecondsValue = ((double) (InputHMSTime % 100)) / 3600.0; 137 /* System.out.printf("calctime : secs=%f\n", SecondsValue); */ 138 double MinutesValue = ((double) ((InputHMSTime / 100) % 100)) / 60.0; 139 /* System.out.printf("calctime : mins=%f\n", MinutesValue); */ 140 double HoursValue = (double) (InputHMSTime / 10000); 141 /* System.out.printf("calctime : hrs=%f\n", HoursValue); */ 142 double PartialDateValue = (HoursValue + MinutesValue + SecondsValue) / 24.0; 143 /* System.out.printf("calctime : pdv=%f\n", PartialDateValue); */ 144 TimeReturnValue = (double) InputJulianDate + (double) PartialDateValue; 145 /* System.out.printf("calctime : retval=%f\n", TimeReturnValue); */ 146 } 147 148 return TimeReturnValue; 149 } 150 151 /** 152 * Convert YYYYMMMDD format character string value to julian date Integer. 153 * 154 * @param InputDateCharString 155 * Character representation of date in YYYYmonDD format where DD 156 * and YYYY are integers and mon is a three character 157 * abbreviation for the month (e.g. 2000MAR13). 158 * 159 * @return Julian Date. Note: value may be {@code 0} if there was a problem. 160 */ 161 public static int cmonth2julian(String InputDateCharString) { 162 /* julian date return value */ 163 int JulianDateReturnValue = 99999; 164 165 int MonthValue = 0; 166 167 String delims = "[A-Z]+"; 168 169 String[] tokens = InputDateCharString.split(delims); 170 171 int YearValue = Integer.parseInt(tokens[0]); 172 int DayValue = Integer.parseInt(tokens[1]); 173 174 Scanner InputValue = new Scanner(InputDateCharString); 175 InputValue.useDelimiter("[^A-Z]+"); 176 String MonthCharString = InputValue.next(); 177 InputValue.close(); 178 179 /* calculate integer month value */ 180 while ((MonthValue < 12) && (!MonthCharString.equals(Months[MonthValue]))) { 181 MonthValue++; 182 } 183 184 /* determine julian date from day/month/year */ 185 if (MonthValue < 12) { 186 JulianDateReturnValue = Functions.idmyyd(DayValue, MonthValue + 1, YearValue); 187 } else { 188 /* Error */ 189 JulianDateReturnValue = 0; 190 } 191 return JulianDateReturnValue; 192 } 193 194 /** 195 * Convert dd/mm/yy to yyyyddd format. 196 * 197 * @param InputDay 198 * Day. 199 * @param InputMonth 200 * Month (integer). 201 * @param InputYear 202 * Year (YY or YYYY). 203 * 204 * @return Julian date or 0 (0=bad input data). 205 */ 206 public static int idmyyd(int InputDay, int InputMonth, int InputYear) { 207 int JulianDate_Return = -1; /* julian date return value */ 208 int DayValue; /* local day value */ 209 210 /* perform a couple quality checks for day/month */ 211 if (((InputDay < 1) || (InputDay > 31)) || ((InputMonth < 1) || (InputMonth > 12))) { 212 JulianDate_Return = 0; 213 } else { 214 DayValue = InputDay + JulianDateMonthArray[InputMonth - 1]; 215 if (InputYear < 1900) { 216 if (InputYear > 70) { 217 InputYear = 1900 + InputYear; 218 } else { 219 InputYear = 2000 + InputYear; 220 } 221 } 222 /* Leap year check */ 223 if (((InputYear % 4) == 0) && (InputMonth > 2)) { 224 DayValue = DayValue + 1; 225 } 226 JulianDate_Return = (InputYear * 1000) + DayValue; 227 } 228 return JulianDate_Return; 229 } 230 231 /** 232 * Convert yyyyddd to dd/mm/yy format. 233 * 234 * @param JulianDateInput 235 * Julian day (yyyyddd). 236 * 237 * @return Array of three values. First value is the day, second value is 238 * the month, and the third value is the year (yyyy). 239 */ 240 public static int[] adt_yddmy(int JulianDateInput) { 241 242 /* local year value */ 243 int YearValue_Local; 244 245 /* local day value */ 246 int DayValue_Local; 247 248 /* local month value */ 249 int MonthValue_Local; 250 251 YearValue_Local = JulianDateInput / 1000; 252 if (YearValue_Local < 1900) { 253 if (YearValue_Local > 70) { 254 YearValue_Local = YearValue_Local + 1900; 255 } else { 256 YearValue_Local = YearValue_Local + 2000; 257 } 258 } 259 DayValue_Local = (JulianDateInput % 1000); 260 for (MonthValue_Local = 0; MonthValue_Local < 13; MonthValue_Local++) { 261 if (DayValue_Local <= JulianDateMonthArray[MonthValue_Local]) { 262 break; 263 } 264 } 265 int Year_Return = YearValue_Local; 266 int Month_Return = MonthValue_Local; 267 int Day_Return = DayValue_Local - JulianDateMonthArray[MonthValue_Local - 1]; 268 if (((YearValue_Local % 4) == 0) && (MonthValue_Local > 2)) { 269 Day_Return = Day_Return + 1; 270 } 271 272 return new int[] { Day_Return, Month_Return, Year_Return }; 273 274 } 275 276 /** 277 * Convert julian date to YYYYMMMDD format for output. 278 * 279 * @param JulianDateInputValue 280 * Julian date representation of date. 281 * 282 * @return Character representation of date in YYYYmonDD format where DD and 283 * YYYY are integers and mon is a three character abbreviation for 284 * the month (e.g. 2000MAR13) 285 */ 286 public static String adt_julian2cmonth(int JulianDateInputValue) { 287 /* calculate date/month/year from julian date */ 288 int ReturnValues[] = Functions.adt_yddmy(JulianDateInputValue); 289 int DayValue = ReturnValues[0]; 290 int MonthValue = ReturnValues[1]; 291 int YearValue = ReturnValues[2]; 292 293 /* form character string representation from various components */ 294 String ReturnDateCharString = String.format("%04d%3s%02d", YearValue, 295 Months[MonthValue - 1], DayValue); 296 297 return ReturnDateCharString; 298 } 299 300 public static double[] distance_angle(double EndLatitudeInput, double EndLongitudeInput, 301 double StartLatitudeInput, double StartLongitudeInput, int UnitFlagID) { 302 double Distance_Intermediate = 0.0; 303 double Distance_Final = 0.0; 304 double Angle_Final = 0.0; 305 double StartLatitude_Radians = StartLatitudeInput * RADIANSCONSTANT; 306 double StartLongitude_Radians = StartLongitudeInput * RADIANSCONSTANT; 307 double EndLatitude_Radians = EndLatitudeInput * RADIANSCONSTANT; 308 double EndLongitude_Radians = EndLongitudeInput * RADIANSCONSTANT; 309 double StartLatitude_Radians_ACOS = Math.cos(StartLatitude_Radians); 310 double StartLongitude_Radians_ACOS = Math.cos(StartLongitude_Radians); 311 double StartLatitude_Radians_ASIN = Math.sin(StartLatitude_Radians); 312 double StartLongitude_Radians_ASIN = Math.sin(StartLongitude_Radians); 313 double EndLatitude_Radians_ACOS = Math.cos(EndLatitude_Radians); 314 double EndLongitude_Radians_ACOS = Math.cos(EndLongitude_Radians); 315 double EndLatitude_Radians_ASIN = Math.sin(EndLatitude_Radians); 316 double EndLongitude_Radians_ASIN = Math.sin(EndLongitude_Radians); 317 double COSCOS_Difference = (EndLatitude_Radians_ACOS * EndLongitude_Radians_ACOS) 318 - (StartLatitude_Radians_ACOS * StartLongitude_Radians_ACOS); 319 double SINCOS_Difference = (EndLatitude_Radians_ACOS * EndLongitude_Radians_ASIN) 320 - (StartLatitude_Radians_ACOS * StartLongitude_Radians_ASIN); 321 double LatitudeSIN_Difference = EndLatitude_Radians_ASIN - StartLatitude_Radians_ASIN; 322 double AdditiveValue = (COSCOS_Difference * COSCOS_Difference) 323 + (SINCOS_Difference * SINCOS_Difference) 324 + (LatitudeSIN_Difference * LatitudeSIN_Difference); 325 Distance_Intermediate = Math.sqrt(AdditiveValue); 326 327 /* Distance_Final is distance in kilometers */ 328 Distance_Final = 2.0 * asin(Distance_Intermediate / 2.0) * EARTHRADIUSKM; 329 330 if (UnitFlagID == 2) { 331 /* Conversion to Miles */ 332 Distance_Final = ((69.0 * Distance_Final) + 55) / 111.0; 333 } 334 if (UnitFlagID == 3) { 335 /* Conversion to Nautical Miles */ 336 Distance_Final = ((60.0 * Distance_Final) + 55) / 111.0; 337 } 338 339 /* Compute Final Angle */ 340 if (Math.abs(Distance_Final) > 0.0001) { 341 Angle_Final = (sin(StartLongitude_Radians - EndLongitude_Radians) * sin((Math.PI / 2.0) 342 - EndLatitude_Radians)) 343 / sin(Distance_Intermediate); 344 } else { 345 Angle_Final = 0.0; 346 } 347 if (Math.abs(Angle_Final) > 1.0) { 348 Angle_Final = a_sign(1.000, Angle_Final); 349 } 350 Angle_Final = Math.asin(Angle_Final) / RADIANSCONSTANT; 351 if (EndLatitude_Radians < StartLatitude_Radians) { 352 Angle_Final = 180.0 - Angle_Final; 353 } 354 if (Angle_Final < 0.0) { 355 Angle_Final = 360.0 + Angle_Final; 356 } 357 358 return new double[] { Distance_Final, Angle_Final }; 359 } 360 361 /** 362 * Calculate a latitude and longitude position from an initial 363 * latitude/longitude and distance/angle values. 364 * 365 * 366 * @param StartLatitudeInput 367 * Initial latitude. 368 * @param StartLongitudeInput 369 * Initial longitude. 370 * @param DistanceInput 371 * Distance from initial position. 372 * @param AngleInput 373 * Angle from initial position. 374 * 375 * Outputs : EndLatitude_Return - derived latitude 376 * EndLongitude_Return - derived longitude 377 */ 378 public static double[] distance_angle2(double StartLatitudeInput, double StartLongitudeInput, 379 double DistanceInput, double AngleInput) { 380 double StartLatitude_Radians = (90.0 - StartLatitudeInput) * RADIANSCONSTANT; 381 double StartLatitude_Radians_Flipped = StartLatitude_Radians; 382 double StartLongitude_Radians = StartLongitudeInput * RADIANSCONSTANT; 383 double ArgumentValue; 384 385 if (StartLatitudeInput < 0.0) { 386 StartLatitude_Radians_Flipped = -(90.0 + StartLatitudeInput) * RADIANSCONSTANT; 387 StartLongitude_Radians = (StartLongitudeInput - 180.0) * RADIANSCONSTANT; 388 AngleInput = 360.0 - AngleInput; 389 } 390 int AngleInput_Integer = (int) AngleInput; 391 double Angle_Radians = -1.0 * ((double) ((540 - AngleInput_Integer) % 360)) 392 * RADIANSCONSTANT; 393 double Distance_Radians = (DistanceInput / 111.1) * RADIANSCONSTANT; 394 double Latitude_AACOS = Math.acos((Math.cos(StartLatitude_Radians) * Math 395 .cos(Distance_Radians)) 396 + (Math.sin(StartLatitude_Radians) * Math.sin(Distance_Radians) * Math 397 .cos(Angle_Radians))); 398 double Longitude_AASIN = 0.0; 399 double ATANValue = 0.0; 400 if (Math.abs(Latitude_AACOS) >= 0.0000001) { 401 ArgumentValue = (Math.sin(Distance_Radians) * Math.sin(Angle_Radians)) 402 / Math.sin(Latitude_AACOS); 403 if (Math.abs(ArgumentValue) > 1.0) { 404 ArgumentValue = a_sign(1.0, ArgumentValue); 405 } 406 Longitude_AASIN = Math.asin(ArgumentValue); 407 ATANValue = (Math.atan(Math.sin(ANGLE90DEGRADIANS - Angle_Radians))) 408 / (Math.tan(ANGLE90DEGRADIANS - Distance_Radians)); 409 if (ATANValue > StartLatitude_Radians_Flipped) { 410 Longitude_AASIN = (2.0 * ANGLE90DEGRADIANS) - Longitude_AASIN; 411 } 412 } 413 Longitude_AASIN = StartLongitude_Radians - Longitude_AASIN; 414 double EndLatitude = 90.0 - (Latitude_AACOS / RADIANSCONSTANT); 415 double EndLongitude = (double) ((int) (10000 * (Longitude_AASIN / RADIANSCONSTANT)) % 3600000) / 10000.0; 416 if (EndLongitude < -180.0) { 417 EndLongitude = EndLongitude + 360.0; 418 } 419 420 return new double[] { EndLatitude, EndLongitude }; 421 422 } 423 424 /** 425 * Determine ocean basin given latitude and longitude position of storm. 426 * Outputs : None Return : basin type (0=atlantic,1=west pacific,2,east 427 * pacific,3=Indian) 428 * 429 * @param LatitudeInput 430 * Latitude. 431 * @param LongitudeInput 432 * Longitude. 433 */ 434 public static int[] adt_oceanbasin(double LatitudeInput, double LongitudeInput) { 435 int BasinID_Local; 436 int DomainID_Local; 437 438 /* flip longitude for McIDAS-V */ 439 LongitudeInput = -1.0 * LongitudeInput; 440 441 if (((LongitudeInput < -180.0) || (LongitudeInput > 180.0)) 442 || ((LatitudeInput < -90.0) || (LatitudeInput > 90.0))) { 443 BasinID_Local = -99; 444 DomainID_Local = -99; 445 } else { 446 if (LatitudeInput >= 0.0) { 447 /* northern hemisphere */ 448 if (LongitudeInput <= -100.0) { 449 BasinID_Local = 1; /* West Pacific */ 450 } else if ((LongitudeInput > -100.0) && (LongitudeInput <= -20.0)) { 451 BasinID_Local = 3; /* Indian */ 452 } else if (LongitudeInput >= 100.0) { 453 BasinID_Local = 2; /* East Pacific */ 454 } else { 455 /* -20 to ~+100 */ 456 if (LatitudeInput > 20.0) { 457 BasinID_Local = 0; /* Atlantic */ 458 } else if (LatitudeInput < 10.0) { 459 if (LongitudeInput < 80.0) { 460 BasinID_Local = 0; /* Atlantic */ 461 } else { 462 BasinID_Local = 2; /* East Pacific */ 463 } 464 } else { 465 /* latitude between 10 and 20 north */ 466 /* 467 * slope of line between (100W,20N) and (80W,10N) is 2 468 * if slope of new point and (100,20) > 2, storm is in 469 * atlantic if slope of new point and (100,20) <= 2, 470 * storm is in pacific 471 */ 472 double AtlEPacDivision = (100.0 - LongitudeInput) / (20.0 - LatitudeInput); 473 if (AtlEPacDivision > 2.0) { 474 BasinID_Local = 0; /* Atlantic */ 475 } else { 476 BasinID_Local = 2; /* East Pacific */ 477 } 478 } 479 } 480 } else { 481 /* southern hemisphere */ 482 if (LongitudeInput <= -135.0) { 483 BasinID_Local = 1; /* West Pacific */ 484 } else if ((LongitudeInput > -135.0) && (LongitudeInput <= -20.0)) { 485 BasinID_Local = 3; /* Indian */ 486 } else if ((LongitudeInput > -20.0) && (LongitudeInput <= 67.0)) { 487 BasinID_Local = 0; /* Atlantic */ 488 } else { 489 BasinID_Local = 2; /* East Pacific */ 490 } 491 } 492 } 493 494 int DomainID_Input = Env.DomainID; 495 if (DomainID_Input == -1) { 496 /* automatically determined storm basin ID */ 497 /* if(BasinID_Local==0) { */ 498 if ((BasinID_Local == 0) || (BasinID_Local == 2)) { 499 DomainID_Local = 0; /* atlantic and East Pac only */ 500 } else { 501 DomainID_Local = 1; /* west pacific and other regions */ 502 } 503 } else { 504 /* manually determined storm basin ID */ 505 DomainID_Local = DomainID_Input; 506 } 507 508 return new int[] { BasinID_Local, DomainID_Local }; 509 510 } 511 512 /** 513 * Calculate slope or y-intercept of all points over SearchTimeInterval 514 * period Inputs : SearchTimeInterval - time period to calculate slope or 515 * y-intercept SlopeInterceptFlag - flag value indicating slope or 516 * y-intercept calculation and parameter to utilize 1 = Final T# (will 517 * return slope) 2 = Adjusted Raw T# (will return slope) 3 = latitude (will 518 * return y-intercept) 4 = longitude (will return y-intercept) 519 * HistoryCurrentPointer_Global contains current analysis information 520 * 521 * @return Slope or Y-Intercept value of line over time period desired 522 */ 523 public static double adt_slopecal(double SearchTimeInterval, int SlopeInterceptFlag) { 524 double Slope; 525 double SlopeYIntReturnValue = 0.0; 526 int CounterMinimum = 4; /* counter minimum value */ 527 int RecDate; 528 int RecTime; 529 int RecLand; 530 int Counter = 0; 531 double XAxisValue = 0.0; 532 double YAxisValue = 0.0; 533 double FirstLonValue = 0.0; 534 double HistoryRecTime; 535 double XMean; 536 double YMean; 537 double XVariance; 538 double YVariance; 539 double CovarianceXY; 540 double RValue; 541 double YIntercept; 542 double SumX = 0.0; /* X-axis sum value */ 543 double SumY = 0.0; /* Y-axis sum value */ 544 double SumSquaresX = 0.0; /* X-axis sum squares value */ 545 double SumSquaresY = 0.0; /* Y-axis sum squares value */ 546 double SumXY = 0.0; /* X*Y sum value */ 547 boolean FoundValidRecordTF = false; /* found valid record logical */ 548 boolean OverWaterTF = false; /* TC over water logical value */ 549 boolean FirstLonValueTF = true; /* 1st lon value for extrap logical */ 550 551 int NumRecsHistory = History.HistoryNumberOfRecords(); 552 553 int ImageDate = History.IRCurrentRecord.date; 554 int ImageTime = History.IRCurrentRecord.time; 555 double CurrentTime = Functions.calctime(ImageDate, ImageTime); 556 double TimeThreshold = CurrentTime - (SearchTimeInterval / 24.0); 557 558 boolean LandFlagTF = Env.LandFlagTF; 559 double InitStrengthValue = Env.InitRawTValue; 560 561 int XInc = 0; 562 while (XInc < NumRecsHistory) { 563 RecDate = History.HistoryFile[XInc].date; 564 RecTime = History.HistoryFile[XInc].time; 565 HistoryRecTime = Functions.calctime(RecDate, RecTime); 566 RecLand = History.HistoryFile[XInc].land; 567 568 if ((HistoryRecTime < CurrentTime) && (HistoryRecTime >= TimeThreshold)) { 569 OverWaterTF = true; 570 if (SlopeInterceptFlag < 3) { 571 if (((LandFlagTF) && (RecLand == 1)) || (InitStrengthValue < 1.0)) { 572 OverWaterTF = false; 573 } 574 } 575 if (OverWaterTF) { 576 XAxisValue = (double) (CurrentTime - HistoryRecTime); 577 if (SlopeInterceptFlag == 1) { 578 YAxisValue = History.HistoryFile[XInc].Tfinal; 579 } 580 if (SlopeInterceptFlag == 2) { 581 YAxisValue = History.HistoryFile[XInc].Traw; 582 } 583 if (SlopeInterceptFlag == 3) { 584 YAxisValue = History.HistoryFile[XInc].latitude; 585 } 586 if (SlopeInterceptFlag == 4) { 587 YAxisValue = History.HistoryFile[XInc].longitude; 588 if (FirstLonValueTF) { 589 FirstLonValue = YAxisValue; 590 FirstLonValueTF = false; 591 } else { 592 if ((FirstLonValue > 100.0) && (YAxisValue < -100.0)) { 593 /* dateline cross W to E */ 594 YAxisValue = YAxisValue + 360.0; 595 } 596 if ((FirstLonValue < -100.0) && (YAxisValue > 100.0)) { 597 /* dateline cross E to W */ 598 YAxisValue = YAxisValue - 360.0; 599 } 600 } 601 } 602 SumX = SumX + XAxisValue; 603 SumY = SumY + YAxisValue; 604 SumSquaresX = SumSquaresX + (XAxisValue * XAxisValue); 605 SumSquaresY = SumSquaresY + (YAxisValue * YAxisValue); 606 SumXY = SumXY + (YAxisValue * XAxisValue); 607 Counter++; 608 FoundValidRecordTF = true; 609 } 610 } else { 611 if (FoundValidRecordTF) { 612 break; 613 } 614 } 615 XInc++; 616 } 617 /* if calculating slope of Final T# values, add in current value */ 618 if (SlopeInterceptFlag <= 2) { 619 CounterMinimum = 6; 620 /* add current record to slope calculation */ 621 if (SlopeInterceptFlag == 1) { 622 YAxisValue = History.IRCurrentRecord.Tfinal; 623 } 624 if (SlopeInterceptFlag == 2) { 625 YAxisValue = History.IRCurrentRecord.Traw; 626 } 627 SumY = SumY + YAxisValue; 628 SumSquaresY = SumSquaresY + (YAxisValue * YAxisValue); 629 Counter++; 630 } 631 632 /* 633 * compute least squares fit of line for data points using the equation 634 * Y = Y* + r(varX/varY)(X - X*) = Mx + B X* = mean of X values (time 635 * values) = XMean Y* = mean of Y values (T# values) = YMean varX = 636 * variance of X = XVariance varY = variance of Y = YVariance r = 637 * covariance - Y*X* = RValue M = slope of line (desired value) = 638 * RValue*(sqrt(YVariance/XVariance)) B = y-intercept = 639 * YMean-(slopecal*XMean) 640 */ 641 /* must have more than CounterMinimum data values to calculate slope */ 642 if (Counter < CounterMinimum) { 643 if ((SlopeInterceptFlag == 3) || (SlopeInterceptFlag == 4)) { 644 Slope = 999.99; 645 } else { 646 Slope = 0.0; 647 } 648 SlopeYIntReturnValue = Slope; 649 } else { 650 XMean = SumX / (double) Counter; 651 YMean = SumY / (double) Counter; 652 XVariance = (SumSquaresX / (double) Counter) - (XMean * XMean); 653 YVariance = (SumSquaresY / (double) Counter) - (YMean * YMean); 654 CovarianceXY = (SumXY / (double) Counter) - (XMean * YMean); 655 RValue = CovarianceXY / Math.sqrt(XVariance * YVariance); 656 if ((Math.abs(XVariance) <= 0.0001) || (Math.abs(YVariance) <= 0.0001)) { 657 Slope = 0.0; 658 } else { 659 Slope = RValue * (Math.sqrt(YVariance / XVariance)); 660 } 661 Slope = (double) ((int) (Slope * 10.0)) / 10.0; 662 if ((SlopeInterceptFlag == 3) || (SlopeInterceptFlag == 4)) { 663 YIntercept = YMean - (Slope * XMean); 664 if (SlopeInterceptFlag == 4) { 665 if (YIntercept < -180.0) 666 YIntercept = YIntercept + 360.0; 667 if (YIntercept > 180.0) 668 YIntercept = YIntercept - 360.0; 669 } 670 /* y-intercept for latitude/longitude extrapolation */ 671 SlopeYIntReturnValue = YIntercept; 672 } else { 673 /* slope for Final T# slope calculation */ 674 SlopeYIntReturnValue = Slope; 675 } 676 } 677 678 return SlopeYIntReturnValue; 679 } 680 681 /** 682 * Obtain pressure or wind value (for Atlantic or West Pacific storms) given 683 * the intensity estimate value. 684 * 685 * @param PressureWindIDValue 686 * Flag for wind (1) or pressure (0) output. 687 * @param CIValue 688 * Current Intensity (CI) value 689 * @param LatitudeInput 690 * Latitude value of analysis 691 * @param LongitudeInput 692 * Longitude value of analysis 693 * 694 * @return Pressure or Wind Speed value. 695 */ 696 public static double adt_getpwval(int PressureWindIDValue, double CIValue, 697 double LatitudeInput, double LongitudeInput) { 698 double PWReturnValue = -999.0; 699 double ROCI_Local = 0.0; 700 double R34_Local = 0.0; 701 int XInc = 2; 702 703 /* use traditional MSLP/Wind-T# conversion */ 704 /* determine correct pressure/wind array bin */ 705 while (((CIValue - 0.001) > PW_TnoValues[XInc]) && (XInc < 82)) { 706 XInc++; 707 } 708 709 int[] ReturnValues = Functions.adt_oceanbasin(LatitudeInput, LongitudeInput); 710 /* int DomainID = Env.DomainID; */ 711 int DomainID = ReturnValues[1]; 712 713 boolean UseCKZTF = Env.UseCKZTF; 714 double CKZGaleRadius34 = Env.CKZGaleRadius; 715 double CKZPenv = Env.CKZPenv; 716 717 /* convert CI value to wind/pressure value */ 718 if (PressureWindIDValue == 1) { 719 PWReturnValue = PW_WindValues[XInc]; /* WIND */ 720 } else { 721 PWReturnValue = PW_PressureValues[DomainID][XInc]; /* PRESSURE */ 722 if (UseCKZTF) { 723 /* 724 * use new Knaff/Zehr T#->MSLP conversion. Need Vmax calculation 725 * and two command line input values 726 */ 727 728 double StormSpeed = 11.0; /* 11 knots -- climo value */ 729 double Vmax_Local = PW_WindValues[XInc]; /* WIND */ 730 731 /* System.out.printf("CKZGaleRadius34=%f\n",CKZGaleRadius34); */ 732 if (CKZGaleRadius34 < 0.0) { 733 ROCI_Local = Math.abs(CKZGaleRadius34); 734 R34_Local = (0.354 * ROCI_Local) + 13.3; 735 } else { 736 R34_Local = CKZGaleRadius34; 737 } 738 739 /* System.out.printf("R34_local=%f\n",R34_Local); */ 740 /* System.out.printf("LatitudeInput=%f\n",LatitudeInput); */ 741 double Vstorm_Local = Vmax_Local - (1.5 * Math.pow(StormSpeed, 0.63)); 742 /* System.out.printf("Vstorm_local=%f\n",Vstorm_Local); */ 743 double EXP_Local = 0.1147 + (0.0055 * Vstorm_Local) 744 - (0.001 * (Math.abs(LatitudeInput) - 25.0)); 745 /* System.out.printf("EXP_local=%f\n",EXP_Local); */ 746 double Vstorm_500 = (R34_Local / 9.0) - 3.0; 747 /* System.out.printf("Vstorm_500=%f\n",Vstorm_500); */ 748 double Vstorm_500c = Vstorm_Local 749 * Math.pow(((66.785 - (0.09102 * Vstorm_Local) + (1.0619 * (Math 750 .abs(LatitudeInput) - 25.0))) / 500.0), EXP_Local); 751 /* System.out.printf("Vstorm_500c=%f\n",Vstorm_500c); */ 752 double S = Math.max((Vstorm_500 / Vstorm_500c), 0.4); 753 /* System.out.printf("S=%f\n",S); */ 754 PWReturnValue = 23.286 - (0.483 * Vstorm_Local) 755 - Math.pow((Vstorm_Local / 24.254), 2) - (12.587 * S) 756 - (0.483 * Math.abs(LatitudeInput)) + CKZPenv; 757 if (PWReturnValue >= (CKZPenv - 1.0)) { 758 PWReturnValue = CKZPenv - 2.0; 759 } 760 /* System.out.printf("PWReturnValue=%f\n",PWReturnValue); */ 761 History.IRCurrentRecord.r34 = (int) CKZGaleRadius34; 762 History.IRCurrentRecord.MSLPenv = (int) CKZPenv; 763 } 764 } 765 766 return PWReturnValue; 767 768 } 769 770 /** 771 * Obtain satellite name given ID number. 772 * 773 * @param SatelliteIDInput 774 * Internal ADT satellite ID number. 775 * 776 * @return Satellite character string ID. 777 */ 778 public static String adt_sattypes(int SatelliteIDInput) { 779 String SatelliteName; 780 781 if (SatelliteIDInput < 0) { 782 SatelliteName = String.format("%s", "MSNG"); 783 } else { 784 SatelliteName = String.format("%s", " OTHER"); 785 if (SatelliteIDInput == 21) 786 SatelliteName = String.format("%s", " GOES1"); 787 if (SatelliteIDInput == 23) 788 SatelliteName = String.format("%s", " GOES2"); 789 if (SatelliteIDInput == 25) 790 SatelliteName = String.format("%s", " GOES3"); 791 if (SatelliteIDInput == 27) 792 SatelliteName = String.format("%s", " GOES4"); 793 if (SatelliteIDInput == 29) 794 SatelliteName = String.format("%s", " GOES5"); 795 if (SatelliteIDInput == 31) 796 SatelliteName = String.format("%s", " GOES6"); 797 if (SatelliteIDInput == 33) 798 SatelliteName = String.format("%s", " GOES7"); 799 if (SatelliteIDInput == 34) 800 SatelliteName = String.format("%s", " FY2B"); 801 if (SatelliteIDInput == 35) 802 SatelliteName = String.format("%s", " FY2C"); 803 if (SatelliteIDInput == 36) 804 SatelliteName = String.format("%s", " FY2D"); 805 if (SatelliteIDInput == 37) 806 SatelliteName = String.format("%s", " FY2E"); 807 if (SatelliteIDInput == 38) 808 SatelliteName = String.format("%s", " FY2F"); 809 if (SatelliteIDInput == 39) 810 SatelliteName = String.format("%s", " FY2G"); 811 if (SatelliteIDInput == 40) 812 SatelliteName = String.format("%s", " FY2H"); 813 if ((SatelliteIDInput >= 42) && (SatelliteIDInput <= 45)) { 814 SatelliteName = String.format("%s", " NOAA"); 815 } 816 if (SatelliteIDInput == 51) 817 SatelliteName = String.format("%s", " MSG1"); 818 if (SatelliteIDInput == 52) 819 SatelliteName = String.format("%s", " MSG2"); 820 if (SatelliteIDInput == 53) 821 SatelliteName = String.format("%s", " MSG3"); 822 if (SatelliteIDInput == 54) 823 SatelliteName = String.format("%s", " MET3"); 824 if (SatelliteIDInput == 55) 825 SatelliteName = String.format("%s", " MET4"); 826 if (SatelliteIDInput == 56) 827 SatelliteName = String.format("%s", " MET5"); 828 if (SatelliteIDInput == 57) 829 SatelliteName = String.format("%s", " MET6"); 830 if (SatelliteIDInput == 58) 831 SatelliteName = String.format("%s", " MET7"); 832 if ((SatelliteIDInput >= 60) && (SatelliteIDInput <= 69)) { 833 SatelliteName = String.format("%s", " NOAA"); 834 } 835 if (SatelliteIDInput == 70) 836 SatelliteName = String.format("%s", " GOES8"); 837 if (SatelliteIDInput == 72) 838 SatelliteName = String.format("%s", " GOES9"); 839 if (SatelliteIDInput == 74) 840 SatelliteName = String.format("%s", "GOES10"); 841 if (SatelliteIDInput == 76) 842 SatelliteName = String.format("%s", "GOES11"); 843 if (SatelliteIDInput == 78) 844 SatelliteName = String.format("%s", "GOES12"); 845 if (SatelliteIDInput == 83) 846 SatelliteName = String.format("%s", " GMS4"); 847 if (SatelliteIDInput == 83) 848 SatelliteName = String.format("%s", " GMS5"); 849 if (SatelliteIDInput == 84) 850 SatelliteName = String.format("%s", "MTSAT1"); 851 if (SatelliteIDInput == 85) 852 SatelliteName = String.format("%s", "MTSAT2"); 853 if ((SatelliteIDInput >= 87) && (SatelliteIDInput <= 95)) { 854 SatelliteName = String.format("%s", " DMSP"); 855 } 856 if (SatelliteIDInput == 95) 857 SatelliteName = String.format("%s", " FY1B"); 858 if (SatelliteIDInput == 96) 859 SatelliteName = String.format("%s", " FY1C"); 860 if (SatelliteIDInput == 97) 861 SatelliteName = String.format("%s", " FY1D"); 862 if ((SatelliteIDInput >= 101) && (SatelliteIDInput <= 171)) { 863 SatelliteName = String.format("%s", " MODIS"); 864 } 865 if (SatelliteIDInput == 180) 866 SatelliteName = String.format("%s", "GOES13"); 867 if (SatelliteIDInput == 182) 868 SatelliteName = String.format("%s", "GOES14"); 869 if (SatelliteIDInput == 184) 870 SatelliteName = String.format("%s", "GOES15"); 871 if (SatelliteIDInput == 186) 872 SatelliteName = String.format("%s", "GOES16"); 873 if ((SatelliteIDInput >= 195) && (SatelliteIDInput <= 196)) { 874 SatelliteName = String.format("%s", " DMSP"); 875 } 876 if (SatelliteIDInput == 230) 877 SatelliteName = String.format("%s", "KLPNA1"); 878 if (SatelliteIDInput == 240) 879 SatelliteName = String.format("%s", "MetOpA"); 880 if (SatelliteIDInput == 241) 881 SatelliteName = String.format("%s", "MetOpB"); 882 if (SatelliteIDInput == 242) 883 SatelliteName = String.format("%s", "MetOpC"); 884 } 885 886 return SatelliteName; 887 888 } 889 890 /** 891 * Derived ATCF file name using various input parameters. 892 * 893 * Inputs : StormName_Input - Storm ID character string value (e.g. 12L) 894 * Outputs : ATCFFileName - Complete file name and path for ATCF file Return 895 * : 0 896 */ 897 public static String adt_atcffilename(String StormNameInput, String StormSiteIDInput) { 898 String ATCFFileName = ""; 899 900 int CurDate = History.IRCurrentRecord.date; 901 int CurTime = History.IRCurrentRecord.time; 902 int ReturnValues[] = Functions.adt_yddmy(CurDate); 903 int DayValue = ReturnValues[0]; 904 int MonthValue = ReturnValues[1]; 905 int YearValue = ReturnValues[2]; 906 907 int HHMMTime = CurTime / 100; 908 ATCFFileName += String.format("%s%5s_%4d%02d%02d%04d_%3s_%3s", StormSiteIDInput, "_DVTO", 909 YearValue, MonthValue, DayValue, HHMMTime, StormNameInput, "FIX"); 910 911 return ATCFFileName; 912 } 913 914 /** 915 * Approximate local zenith angle (angle from vertical to the satellite). 916 * 917 * @param rlat 918 * Latitude of target position on earth. 919 * @param rlon 920 * Longitude of target position on earth. 921 * @param plat 922 * Latitude of sub-satellite point. 923 * @param plon 924 * Longitude of sub-satellite point. 925 * 926 * @return local zenith angle 927 */ 928 public static double adt_XLZA(double rlat, double rlon, double plat, double plon) { 929 /* mean earth radius */ 930 double r = 6371.229; 931 /* mean distance from satellite to earth center */ 932 double sd = r + 35790.0; 933 934 rlat = Math.abs(rlat); 935 /* For the longitude in East of Greenwich */ 936 if (rlon < 0) 937 rlon = rlon + 360.0; 938 /* rlon = Math.abs(rlon); is this even needed with above + 360.0? */ 939 double dif = Math.abs(plon - rlon); 940 941 double cosal = Math.cos(deg_to_rad(rlat)) * Math.cos(deg_to_rad(dif)); 942 double dp = Math.sqrt(Math.pow(r, 2) + Math.pow(sd, 2) - (2.0 * r * sd * cosal)); 943 double yy1 = Math.sqrt(1 - Math.pow(cosal, 2)) / cosal; 944 double alpha = Math.atan(yy1); 945 double zz = sd * Math.sin(alpha) / dp; 946 double yy2 = zz / (Math.sqrt(1.0 - Math.pow(zz, 2))); 947 double zaz = rad_to_deg(Math.atan(yy2)); 948 949 return zaz; 950 951 } 952 953}