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.data.cyclone; 030 031import java.util.ArrayList; 032import java.util.Iterator; 033import java.util.List; 034 035/** 036 * Created by IntelliJ IDEA. User: yuanho Date: Feb 26, 2009 Time: 10:02:45 AM 037 * To change this template use File | Settings | File Templates. 038 */ 039 040public class StormAODTSceneType { 041 042 /** _more_ */ 043 static int maxsec = 24; 044 045 /** _more_ */ 046 static int maxsecA = 2000; 047 048 /** _more_ */ 049 static double PI = 3.14159265358979323846; 050 051 /** 052 * _more_ 053 * 054 * @param odtcurrent 055 * _more_ 056 * @param areadata 057 * _more_ 058 * 059 * @return _more_ 060 */ 061 static float[] aodtv72_calcscene(StormAODTInfo.IRData odtcurrent, 062 StormAODTInfo.DataGrid areadata) 063 /* 064 * Perform Fast Fourier Transform (FFT) analysis and determine scene type 065 * via empirically defined threshold values. Inputs : global structure 066 * odtcurrent_v72 containing current intensity values Outputs : various 067 * elements within structure odtcurrent_v72 Return : -41 : error with FFT 068 * routine -51 : cloud temperature value <-100C or >+40C 0 : o.k. 069 */ 070 { 071 072 int nbin = 64, iok, a, b, bptr; 073 int ixx, iyy, izz, iscene, idxcnt, maxsecd2, sxscnt; 074 int[] sectorcnt; 075 int cnteye; 076 float[] bd, cbd, tbd; 077 078 float radb, rade, radeye, teye, dx2, eyecnt = 0, rngcnt = 0; 079 float xangle, xdist, xtemp, slice; 080 float[] sectordiffa, sectormin; 081 float[][] sector; 082 float[] eyearr, sxs; 083 float sectang1, sectang2; 084 float[] avex, stdvx, skewx, xs, i2; 085 float Aaveext, Astdvext, Askewext, Aavesym, Astdvsym, Askewsym; 086 float Eaveeye, Estdveye, Eskeweye; 087 float alsdist, alsradb, alsrade, alssum, alst, alscnt; 088 StormAODTInfo.RingData tcirc; 089 090 float eyetemp; 091 float eyecdosize, rmw; 092 int eyecat; 093 094 /* initialize temperature histogram bin values */ 095 bd = new float[534]; 096 cbd = new float[534]; 097 tbd = new float[534]; 098 for (iyy = 0; iyy < nbin; iyy++) { 099 bd[iyy] = 26.0f - (float) iyy * 2.0f; 100 } 101 /* 102 * set up arrays for FFT ananlysis. iscene=0 will perform FFT for cloud 103 * top region while iscene=1 will perform FFT for eye region 104 */ 105 int cursorx = areadata.numx / 2; 106 int cursory = areadata.numy / 2; 107 List<StormAODTInfo.RingData> tcircfirst = aodtv72_readcirc(cursorx, 108 cursory, areadata); 109 110 for (iscene = 0; iscene <= 1; iscene++) { 111 for (iyy = 0; iyy < nbin; iyy++) { 112 cbd[iyy] = 0.0f; 113 tbd[iyy] = 0.0f; 114 } 115 116 /* define start and end radii values */ 117 if (iscene == 0) { 118 /* CLOUD TOP */ 119 radb = (float) StormAODTInfo.kstart_v72; 120 rade = (float) StormAODTInfo.kend_v72; 121 } else { 122 /* EYE REGION */ 123 radb = 0.0f; 124 rade = (float) StormAODTInfo.kstart_v72; 125 } 126 127 /* load arrays for FFT analysis */ 128 Iterator iit = tcircfirst.iterator(); 129 130 while (iit.hasNext()) { 131 tcirc = (StormAODTInfo.RingData) iit.next(); 132 if ((tcirc.dist >= radb) && (tcirc.dist <= rade)) { 133 teye = (tcirc.temp - 273.16f); 134 for (ixx = 0; ixx < (nbin - 1); ixx++) { 135 if ((teye <= bd[ixx]) && (teye >= bd[ixx + 1])) { 136 cbd[ixx] = cbd[ixx] + 1.0f; 137 tbd[ixx] = tbd[ixx] + teye; 138 } 139 } 140 } 141 } 142 143 /* perform FFT analysis */ 144 float[] fftOut = aodtv72_fft(cbd); 145 if (fftOut == null) { 146 return null; // return -41; 147 } 148 149 /* assign variables based upon region being analyzed */ 150 if (iscene == 0) { 151 rngcnt = fftOut[1]; /* CLOUD TOP */ 152 } else { 153 eyecnt = fftOut[1]; /* EYE REGION */ 154 } 155 } 156 157 /* 158 * compute various cloud and eye region parameters for classification 159 * scheme 160 */ 161 162 /* allocate memory for arrays */ 163 164 eyearr = new float[maxsecA]; 165 sectorcnt = new int[maxsec]; 166 sectormin = new float[maxsec]; 167 sector = new float[maxsec][maxsecA]; 168 for (ixx = 0; ixx < maxsec; ixx++) { 169 sectorcnt[ixx] = 0; 170 sectormin[ixx] = 999.0f; 171 for (iyy = 0; iyy < maxsecA; iyy++) { 172 sector[ixx][iyy] = -999.99f; 173 } 174 } 175 176 /* load array for analysis */ 177 radb = (float) StormAODTInfo.kstart_v72; 178 rade = (float) StormAODTInfo.kend_v72; 179 radeye = (float) StormAODTInfo.kstart_v72; 180 Iterator it = tcircfirst.iterator(); 181 182 slice = 360.0f / (float) maxsec; 183 cnteye = 0; 184 while (it.hasNext()) { 185 tcirc = (StormAODTInfo.RingData) it.next(); 186 xangle = tcirc.angle; 187 xdist = tcirc.dist; 188 xtemp = tcirc.temp - 273.16f; 189 if (xangle == 360.0) { 190 xangle = 0.0f; 191 } 192 ixx = 0; 193 if ((xdist >= radb) && (xdist <= rade)) { 194 while (ixx < maxsec) { 195 sectang1 = (float) Math.max(0.0, (float) ixx * slice); 196 sectang2 = (float) Math.min(360.0, (float) (ixx + 1) 197 * slice); 198 if ((xangle >= sectang1) && (xangle < sectang2)) { 199 sector[ixx][sectorcnt[ixx]] = xtemp; 200 sectorcnt[ixx]++; 201 ixx = maxsec; 202 } else { 203 ixx++; 204 } 205 } 206 } 207 if ((xdist >= 0) && (xdist <= radeye)) { 208 eyearr[cnteye] = xtemp; 209 cnteye++; 210 } 211 /* 212 * the following will count all values w/in the different BD curve 213 * ranges for the whole cloud region so we can examine the structure 214 * of the cloud region later 215 */ 216 // tcirc=tcirc->nextrec; 217 } 218 219 /* 220 * position annulus at CW max temp distance and determine mean temp w/in 221 * +/- 40km from this distance. If dist is less than 68km from center, 222 * annulus will start at 28km 223 */ 224 alscnt = 0; 225 alssum = 0.0f; 226 alsdist = (float) odtcurrent.cwring; 227 alsradb = Math.max(28.0f, alsdist - 40.0f); 228 alsrade = Math.max(108.0f, alsdist + 40.0f); 229 it = tcircfirst.iterator(); 230 while (it.hasNext()) { 231 tcirc = (StormAODTInfo.RingData) it.next(); 232 xdist = tcirc.dist; 233 xtemp = tcirc.temp - 273.16f; 234 if ((xdist >= alsradb) && (xdist <= alsrade)) { 235 alssum = alssum + xtemp; 236 alscnt++; 237 } 238 // tcirc=tcirc->nextrec; 239 } 240 alst = alssum / (float) alscnt; 241 // odtcurrent.cloudt=alst; 242 // if ((odtcurrent.cloudt < -100.0) || (odtcurrent.cloudt > 40.0)) { 243 if ((alst < -100.0) || (alst > 40.0)) { 244 return null; // return -51; 245 } 246 247 /* 248 * determine "minimum" temperature for each sector. This is not the 249 * actual lowest temperature, but instead is the temperature value of 250 * the coldest 90% of the values within the sector. This will, 251 * hopefully, eliminate some outliers 252 */ 253 254 /* allocate memory */ 255 xs = new float[maxsecA]; 256 i2 = new float[maxsecA]; 257 sxs = new float[maxsecA]; 258 259 for (ixx = 0; ixx < maxsec; ixx++) { 260 sxscnt = sectorcnt[ixx]; 261 for (iyy = 0; iyy < sxscnt; iyy++) { 262 sxs[iyy] = sector[ixx][iyy]; 263 } 264 aodtv72_xxsort(sxs, xs, i2, sxscnt); 265 izz = (int) (sxscnt - (sxscnt * .1f)); 266 sectormin[ixx] = xs[izz]; 267 } 268 269 /* free memory */ 270 271 /* determine averages, standard deviations and skews for each sector */ 272 /* allocate memory */ 273 avex = new float[maxsec]; 274 stdvx = new float[maxsec]; 275 skewx = new float[maxsec]; 276 sectordiffa = new float[maxsec]; 277 278 for (ixx = 0; ixx < maxsec; ixx++) { 279 float[] out = aodtv72_calcskew(sector[ixx], sectorcnt[ixx]); 280 avex[ixx] = out[0]; 281 stdvx[ixx] = out[1]; 282 skewx[ixx] = out[2]; 283 } 284 285 float[] out = aodtv72_calcskew(avex, maxsec); 286 Aaveext = out[0]; 287 Astdvext = out[1]; 288 Askewext = out[2]; 289 // odtcurrent.cloudt2=Aaveext; 290 291 /* these are used to examine symmetry of convection */ 292 maxsecd2 = maxsec / 2; 293 for (ixx = 0; ixx < maxsecd2; ixx++) { 294 sectordiffa[ixx] = Math.abs(avex[ixx] - avex[maxsecd2 + ixx]); 295 } 296 out = aodtv72_calcskew(sectordiffa, maxsecd2); 297 Aavesym = out[0]; 298 Astdvsym = out[1]; 299 Askewsym = out[2]; 300 /* these are used to examine properties of the eye region */ 301 out = aodtv72_calcskew(eyearr, cnteye); 302 Eaveeye = out[0]; 303 Estdveye = out[1]; 304 Eskeweye = out[2]; 305 306 // odtcurrent.eyestdv=Estdveye; 307 // odtcurrent.cloudsymave=Aavesym; 308 // odtcurrent.eyefft=(int)eyecnt; 309 // odtcurrent.cloudfft=(int)rngcnt; 310 311 /* free memory */ 312 313 /* assign scenetype value in structure odtcurrent_v72 */ 314 315 return new float[] { alst, Aaveext, Estdveye, Aavesym, eyecnt, rngcnt }; 316 317 // return odt; 318 } 319 320 /** 321 * _more_ 322 * 323 * @param odtcurrent 324 * _more_ 325 * @param rmwsizeman 326 * _more_ 327 * @param areadata_v72 328 * _more_ 329 * @param osstr_v72 330 * _more_ 331 * @param osearch_v72 332 * _more_ 333 * 334 * @return _more_ 335 */ 336 337 static float[] aodtv72_classify(StormAODTInfo.IRData odtcurrent, 338 int rmwsizeman, StormAODTInfo.DataGrid areadata_v72, 339 float osstr_v72, boolean osearch_v72) 340 341 /* 342 * Classify scene type based on FFT analysis and histogram temperatures 343 * using empirically defined threshold values. Inputs : global structure 344 * odtcurrent_v72 containing current image analysis Outputs : scenetype in 345 * structure odtcurrent_v72 is modified 346 * 347 * SCENE TYPES : EYE REGION CLOUD REGION 0 - clear 0 - uniform 1 - pinhole 1 348 * - embedded center 2 - large 2 - irregular cdo 3 - none 3 - curved band 4 349 * - shear 350 */ 351 { 352 353 int iok, ixx, iyy, cloudfft = 0, cloudcat = 0, eyefft, eyecat = 0, cwtcat = 0, diffcat; 354 int cbring, cbringval, diffcloudcat; 355 int sceneeye, scenecloud, spiral = 0, spiralmax; 356 int lastcscene, lastescene, cbringvalmax; 357 float xlat, xlon, tempval, lasttno, lasttno12; 358 float eyetemp, cloudtemp, cloudcwt, eyestdv, cloudsyma; 359 float essizeDG = 0.0f, esdiffDG = 0.0f; 360 float lat1, lon1, lat2, lon2; 361 float eyecdosize, diffcloudt, rmw; 362 float cbringlatmax, cbringlonmax; 363 double curtime, curtimem12, xtime, lastvalidtime; 364 365 boolean cbfound, cbscene, cbgray, shear, irrcdo, landcheck; 366 float[] cdodiam = { 0.0f, 85.0f, 140.0f, 195.0f, 250.0f, 999.0f }; 367 float xpart, ddvor; 368 int cdocat; 369 370 boolean embdd = false, checkembc = false; 371 float Ea, Eb, Ec, Ed, Ee, eyeval, Ca, Cb, Cc, Cd, Ce, cloudval; 372 float diffeyecloudt, eyecloudcatdiff, eyecwtcatdiff, cloudcatdiff; 373 float eyepart = 0, cloudpart = 0, cwtpart = 0, diffcat2, lastvalidtno, lastr9, maxtno; 374 float[] cdomax = { 0.0f, 0.4f, 0.4f, 0.5f, 0.5f, 0.6f, 0.6f }; 375 int diffeyecloudcat; 376 boolean foundeye, foundm12; 377 378 eyetemp = odtcurrent.eyet; 379 eyefft = odtcurrent.eyefft; 380 eyestdv = odtcurrent.eyestdv; 381 cloudtemp = odtcurrent.cloudt; 382 cloudcwt = odtcurrent.cwcloudt; 383 cloudfft = odtcurrent.cloudfft; 384 cloudsyma = odtcurrent.cloudsymave; 385 xlat = odtcurrent.latitude; 386 xlon = odtcurrent.longitude; 387 388 // odtcurrent.eyestdv=Estdveye; 389 // odtcurrent.cloudsymave=Aavesym; 390 // odtcurrent.eyefft=(int)eyecnt; 391 // odtcurrent.cloudfft=(int)rngcnt; 392 393 for (ixx = 0; ixx < 10; ixx++) { 394 /* compute cloud category */ 395 if ((cloudtemp <= StormAODTInfo.ebd_v72[ixx]) 396 && (cloudtemp > StormAODTInfo.ebd_v72[ixx + 1])) { 397 cloudcat = ixx; 398 xpart = (float) ((cloudtemp - StormAODTInfo.ebd_v72[cloudcat]) / (StormAODTInfo.ebd_v72[cloudcat + 1] - StormAODTInfo.ebd_v72[cloudcat])); 399 if (cloudcat == 0) { 400 xpart = 0.0f; 401 } 402 cloudpart = cloudcat + xpart; 403 } 404 /* compute eye category */ 405 if ((eyetemp <= StormAODTInfo.ebd_v72[ixx]) 406 && (eyetemp > StormAODTInfo.ebd_v72[ixx + 1])) { 407 eyecat = ixx; 408 xpart = (float) ((eyetemp - StormAODTInfo.ebd_v72[eyecat]) / (StormAODTInfo.ebd_v72[eyecat + 1] - StormAODTInfo.ebd_v72[eyecat])); 409 if (eyecat == 0) { 410 xpart = 0.0f; 411 } 412 eyepart = eyecat + xpart; 413 } 414 /* compute C-W eye category */ 415 if ((cloudcwt <= StormAODTInfo.ebd_v72[ixx]) 416 && (cloudcwt > StormAODTInfo.ebd_v72[ixx + 1])) { 417 cwtcat = ixx; 418 xpart = (float) ((cloudcwt - StormAODTInfo.ebd_v72[cwtcat]) / (StormAODTInfo.ebd_v72[cwtcat + 1] - StormAODTInfo.ebd_v72[cwtcat])); 419 if (cwtcat == 0) { 420 xpart = 0.0f; 421 } 422 cwtpart = cwtcat + xpart; 423 } 424 } 425 426 /* printf("EYE = temp=%f cat=%d part=%f \n",eyetemp,eyecat,eyepart); */ 427 /* 428 * printf("CLD = temp=%f cat=%d part=%f \n",cloudtemp,cloudcat,cloudpart) 429 * ; 430 */ 431 /* printf("CWT = temp=%f cat=%d part=%f \n",cloudcwt,cwtcat,cwtpart); */ 432 diffcat = Math.max(0, Math.max(cloudcat, cwtcat) - eyecat); 433 diffcloudt = cloudtemp - cloudcwt; 434 diffeyecloudt = eyetemp - cloudtemp; 435 eyecwtcatdiff = cwtpart - eyepart; 436 eyecloudcatdiff = cloudpart - eyepart; 437 cloudcatdiff = cloudpart - cwtpart; 438 diffcloudcat = cloudcat - cwtcat; 439 diffeyecloudcat = cloudcat - eyecat; 440 diffcat2 = eyetemp - (Math.min(cloudtemp, cloudcwt)); 441 442 curtime = odtcurrent.date; // aodtv72_calctime(odtcurrent.date, 443 // odtcurrent.time); 444 // curtimem12 = curtime - 0.5; 445 446 /* determine last Final T# value for curved band/other scene check */ 447 foundeye = false; 448 maxtno = 0.0f; 449 lastr9 = 0; 450 lasttno = maxtno; 451 lastvalidtno = lasttno; 452 if (true) { // (odthistoryfirst_v72==0)||(strlen(hfile_v72)==0)) { 453 foundm12 = true; 454 lastescene = 3; 455 lastr9 = 1; 456 if ((cwtpart < 3.5) && (osstr_v72 < 3.5)) { 457 lastcscene = 3; 458 lasttno12 = osstr_v72; 459 } else { 460 lastcscene = 0; 461 lasttno12 = Math.max(osstr_v72, 4.0f); 462 } 463 } /* 464 * else { odthistory=odthistoryfirst_v72; foundm12=FALSE; lastcscene=3; 465 * while(odthistory!=0) { 466 * xtime=aodtv72_calctime(odthistory->IR.date,odthistory->IR.time); 467 * landcheck=true; 468 * if(((oland_v72)&&(odthistory->IR.land==1))||(odthistory 469 * ->IR.Traw<1.0)) landcheck=FALSE; if((xtime<curtime)&&(landcheck)) { 470 * lastvalidtime=xtime; if((xtime>=curtimem12)&&(!foundm12)) { 471 * lasttno12=odthistory->IR.Tfinal; foundm12=true; } 472 * lasttno=odthistory->IR.Tfinal; lastcscene=odthistory->IR.cloudscene; 473 * lastescene=odthistory->IR.eyescene; if(lastescene<=2) foundeye=true; 474 * if((lastcscene==4)&&(lastescene==3)) foundeye=FALSE; 475 * lastvalidtno=lasttno; lastr9=odthistory->IR.rule9; if(lasttno>maxtno) 476 * maxtno=lasttno; } else { if(!landcheck) { 477 * 478 * if((xtime-lastvalidtime)>0.5) { foundeye=FALSE; 479 * lasttno=lastvalidtno-(1.0*(xtime-lastvalidtime)); 480 * 481 * } } } odthistory=odthistory->nextrec; } 482 * 483 * if(!foundm12) lasttno12=lasttno; } 484 */ 485 /* printf("foundm12=%d\n",foundm12); */ 486 /* printf("lasttno12=%f\n",lasttno12); */ 487 488 /* NEW SCENE IDENTIFICATION SCHEME */ 489 /* NEW EYE SCENE THRESHOLD DETERMINATION */ 490 Ec = 0.0f; 491 Ee = 0.0f; 492 Ea = 1.0f - ((eyefft - 2) * 0.1f); 493 Eb = -(eyepart * 0.5f); 494 if (eyestdv > 10.0) { 495 Ec = 0.50f; 496 } 497 /* Ed=eyecloudcatdiff*0.75; */ 498 Ed = (eyecloudcatdiff * 0.25f) + (eyecwtcatdiff * 0.50f); /* new */ 499 /* printf("Ed=%f\n",Ed); */ 500 /* printf("maxtno=%f Ee=%f\n",maxtno,Ee); */ 501 if ((foundm12) && (lastescene < 3) && (maxtno > 5.0)) { 502 Ee = Ee + 0.25f; 503 } 504 /* printf("Ee=%f\n",Ee); */ 505 /* if(lasttno12<3.5) Ee=Ee-1.0; */ 506 if (lasttno12 <= 4.5) { 507 Ee = Math.max(-1.0f, lasttno12 - 4.5f); /* new */ 508 } 509 /* printf("lasttno12=%f Ee=%f\n",lasttno12,Ee); */ 510 if ((lastr9 > 0) && (lasttno < 4.0)) { 511 Ee = Ee - 0.5f; /* new */ 512 } 513 /* printf("Ee=%f\n",Ee); */ 514 eyeval = Ea + Eb + Ec + Ed + Ee; 515 sceneeye = 3; /* NO EYE */ 516 if (eyeval >= 0.50) { 517 sceneeye = 0; /* EYE */ 518 } 519 /* 520 * if(eyeval>=0.00) sceneeye=5; / OBSCURED EYE / if(eyeval>=1.50) 521 * sceneeye=4; / RAGGED EYE / if(eyeval>=3.50) sceneeye=0; / CLEAR EYE / 522 */ 523 524 /* printf("eyeval= %f sceneeye=%d \n",eyeval,sceneeye); */ 525 float odtcurrent_v72IRRMW; 526 eyecdosize = 0.0f; 527 if (rmwsizeman > 0) { 528 /* printf("RMW SIZE=%d\n",rmwsizeman_v72); */ 529 odtcurrent_v72IRRMW = (float) rmwsizeman; 530 eyecdosize = (float) rmwsizeman - 1.0f; /* manually input eye size */ 531 } else { 532 float[] out = aodtv72_rmw(odtcurrent, areadata_v72); 533 rmw = out[0]; 534 eyecdosize = out[1]; 535 odtcurrent_v72IRRMW = rmw; 536 } 537 538 /* LARGE EYE CHECKS */ 539 if ((sceneeye == 0) && (eyecdosize >= 45.0)) { 540 sceneeye = 2; /* large eye */ 541 } 542 543 /* NEW CLOUD SCENE THRESHOLD DETERMINATION */ 544 shear = false; 545 irrcdo = false; 546 cbscene = true; 547 cbgray = true; 548 549 Cc = 0.0f; 550 Cd = 0.5f; /* changed to 0.5 */ 551 Ce = 0.0f; 552 Ca = cwtpart * 0.25f; 553 Cb = cloudpart * 0.25f; 554 if (cloudfft <= 2) { 555 Cc = Math.min(1.50f, cwtpart * 0.25f); 556 } 557 if (lastcscene >= 3) { 558 Cd = -0.50f; 559 } 560 /* printf("cwtpart=%f lasttno12=%f\n",cwtpart,lasttno12); */ 561 if (cwtpart > 2.0) { /* new */ 562 if (lasttno12 >= 2.5) { 563 if (sceneeye == 0) { 564 Ce = Math.min(1.00f, lasttno12 - 2.5f); 565 } 566 if (lasttno12 >= 3.5) { 567 Ce = Ce + 1.00f; 568 } 569 } 570 if ((foundm12) && (foundeye)) { 571 Ce = Ce + 1.25f; 572 } 573 } 574 cloudval = Ca + Cb + Cc + Cd + Ce; 575 if (cloudval < 0.0) { 576 shear = true; /* SHEAR */ 577 } 578 if (cloudval >= 0.00) { 579 cbscene = true; /* CURVED BAND (gray) */ 580 } 581 if (cloudval >= 1.00) { 582 cbscene = true; /* CURVED BAND (gray) */ 583 /* check for irregular CDO */ 584 if ((diffcat2 < 0.0) && (cloudsyma > 40.0)) { 585 irrcdo = true; /* IRREGULAR CDO */ 586 } 587 } 588 if ((cloudval >= 2.00) && (cloudval < 3.00)) { 589 cbscene = true; /* CURVED BAND (gray) */ 590 /* check for irregular CDO */ 591 /* if((diffcloudcat<0)&&(diffcloudt>8.0)&&(cloudsyma>30.0)) { */ 592 if ((diffcat2 < 0.0) && (cloudsyma > 30.0)) { 593 irrcdo = true; /* IRREGULAR CDO */ 594 } 595 if (cwtcat >= 3) { 596 /* if xcwt>3.0 try black/white CB check */ 597 if ((diffcloudcat > 0) && (diffcloudt < -8.0)) { 598 cbgray = false; /* CURVED BAND (black/white) */ 599 } 600 /* check for large/ragged eye */ 601 if ((sceneeye == 0) 602 || ((eyepart > 1.00) && (diffeyecloudcat >= 2.00))) { 603 cbscene = false; /* EYE */ 604 } 605 /* check for CDO */ 606 if ((cloudcatdiff <= 0.0) && (eyecwtcatdiff < 1.00)) { 607 cbscene = false; /* CDO */ 608 } 609 } 610 } 611 if (cloudval >= 3.00) { 612 cbscene = false; /* CDO */ 613 /* check for irregular CDO */ 614 if ((diffcloudcat < 0) && (diffcloudt > 8.0) && (cloudsyma > 30.0)) { 615 irrcdo = true; /* IRREGULAR CDO */ 616 cbscene = true; 617 } 618 } 619 /* EMBEDDED CENTER CHECK */ 620 if ((cloudtemp < cloudcwt) && (cloudcwt < eyetemp)) { 621 checkembc = true; 622 } 623 if ((!cbscene) && (checkembc)) { 624 tempval = (float) StormAODTInfo.ebd_v72[cwtcat + 1] + 273.16f; 625 float[] out = aodtv72_logspiral(xlat, xlon, tempval, 1, odtcurrent, 626 areadata_v72); 627 spiral = (int) out[0]; 628 if ((spiral >= 8) && (spiral < 20)) { 629 embdd = true; 630 } 631 /* printf(" EMBDD : cwtcat=%d spiral=%d \n",cwtcat,spiral); */ 632 } 633 634 /* 635 * printf("cloudval= %f shear=%d cbscene=%d cbgray=%d irrcdo=%d \n",cloudval 636 * ,shear,cbscene,cbgray,irrcdo); 637 */ 638 // (void)aodtv72_julian2cmonth(odtcurrent_v72IR.date,cdate); 639 /* 640 * printf("%9s %6d %4.1f %4.1f %2d %2d %5.1f %5.1f %5.2f %5.2f %5.2f 641 * %5.2f %5.2f %5.2f %5.2f %2d %2d %4.1f %4.1f %5.2f %5.2f %5.2f %5.2f 642 * %5.2f %5.2f %6.2f %7.2f %3.1f \n",cdate,odtcurrent_v72->IR.time, 643 * cloudpart 644 * ,cwtpart,cloudfft,diffcloudcat,diffcloudt,cloudsyma,Ca,Cb,0.0 645 * ,Cc,Cd,Ce,cloudval, 646 * eyefft,diffeyecloudcat,eyepart,eyestdv,Ea,Eb,Ec,Ed 647 * ,Ee,eyeval,xlat,xlon,lasttno); 648 */ 649 650 /* CLASSIFY CLOUD REGION */ 651 cbring = 0; 652 cbringval = 0; 653 cbringvalmax = 0; 654 cbringlatmax = xlat; 655 cbringlonmax = xlon; 656 657 // L100: 658 if (cbscene) { 659 if (shear) { 660 sceneeye = 3; /* NO EYE */ 661 scenecloud = 4; /* SHEAR */ 662 tempval = (float) StormAODTInfo.ebd_v72[3] + 273.16f; 663 float[] out = aodtv72_cdoshearcalc(xlat, xlon, tempval, 3, 664 areadata_v72); 665 eyecdosize = Math.max(4.0f, out[0]); 666 } else if (irrcdo) { 667 sceneeye = 3; /* NO EYE */ 668 scenecloud = 2; /* IRREGULAR CDO */ 669 } else { 670 cbfound = false; 671 boolean cont = true; 672 // L200: 673 if (cbgray) { 674 /* perform Curved Band analysis */ 675 ixx = 4; /* start with LIGHT GRAY */ 676 while ((ixx >= 2) && (!cbfound) && cont) { 677 tempval = (float) StormAODTInfo.ebd_v72[ixx] + 273.16f; 678 float[] out0 = aodtv72_logspiral(xlat, xlon, tempval, 679 1, odtcurrent, areadata_v72); 680 spiral = (int) out0[0]; 681 /* printf("BD level=%d spiral=%d\n",ixx,spiral); */ 682 if ((out0[0] >= 8) || (ixx == 2)) { /* 683 * 10 = .375% -- 10 684 * ==> 9 arcs of 15 685 * degrees 686 */ 687 if (spiral > 25) { 688 if (ixx == 4) { 689 cbgray = false; 690 // goto L200; 691 cbscene = false; 692 ixx = 6; 693 while ((ixx > 4) && (!cbfound)) { 694 tempval = (float) StormAODTInfo.ebd_v72[ixx] + 273.16f; 695 float[] out = aodtv72_logspiral(xlat, 696 xlon, tempval, 1, odtcurrent, 697 areadata_v72); 698 spiral = (int) out[0]; 699 if ((spiral >= 9) && (spiral <= 25)) { 700 cbfound = true; 701 } else { 702 ixx--; 703 } 704 } 705 } else { 706 ixx = 0; 707 } 708 } else { 709 if ((ixx == 2) && (spiral < 7)) { /* 710 * 7 = .25% -- 711 * 10 ==> 6 arcs 712 * of 15 degrees 713 */ 714 /* probably shear */ 715 cbfound = false; 716 shear = true; 717 // goto L100; 718 sceneeye = 3; /* NO EYE */ 719 scenecloud = 4; /* SHEAR */ 720 tempval = (float) StormAODTInfo.ebd_v72[3] + 273.16f; 721 float[] out = aodtv72_cdoshearcalc(xlat, 722 xlon, tempval, 3, areadata_v72); 723 eyecdosize = Math.max(4.0f, out[0]); 724 cont = false; 725 } else { 726 cbfound = true; 727 } 728 } 729 } else { 730 ixx--; 731 } 732 } 733 } else { 734 /* try BLACK and WHITE rings */ 735 cbscene = false; 736 ixx = 6; 737 while ((ixx > 4) && (!cbfound)) { 738 tempval = (float) StormAODTInfo.ebd_v72[ixx] + 273.16f; 739 float[] out = aodtv72_logspiral(xlat, xlon, tempval, 1, 740 odtcurrent, areadata_v72); 741 spiral = (int) out[0]; 742 if ((spiral >= 9) && (spiral <= 25)) { 743 cbfound = true; 744 } else { 745 ixx--; 746 } 747 } 748 } 749 if (cbfound && cont) { 750 /* found curved band scenes */ 751 cbring = ixx; 752 cbringval = spiral; 753 sceneeye = 3; /* NO EYE */ 754 scenecloud = 3; /* CURVED BAND */ 755 /* 756 * search for maximum curved band analysis location within 757 * 1-degree box 758 */ 759 tempval = (float) StormAODTInfo.ebd_v72[cbring] + 273.16f; 760 if (osearch_v72) { 761 float[] out = aodtv72_logspiral(xlat, xlon, tempval, 2, 762 odtcurrent, areadata_v72); 763 spiralmax = (int) out[0]; 764 lat2 = out[1]; 765 lon2 = out[2]; 766 767 cbringvalmax = (int) out[0]; 768 cbringlatmax = out[1]; 769 cbringlonmax = out[2]; 770 } 771 } else if (cont) { 772 /* 773 * did not find curved band scenes, mark as non-eye/eye 774 * scene 775 */ 776 scenecloud = 0; 777 cbscene = false; 778 embdd = false; /* redundant declaration */ 779 // goto L100; 780 scenecloud = 0; /* UNIFORM */ 781 if (embdd) { 782 scenecloud = 1; /* EMBEDDED CENTER */ 783 } 784 /* PINHOLE EYE TEST */ 785 /* 786 * printf("sceneeye=%d diffeyecloudcat=%d eyefft=%d cwtpart=%f scenecloud=%d cloudfft=%d lasttno12=%f\n" 787 * , 788 * sceneeye,diffeyecloudcat,eyefft,cwtpart,scenecloud,cloudfft 789 * ,lasttno12); 790 */ 791 if ((rmwsizeman > 0) && (rmwsizeman < 5.0)) { 792 sceneeye = 1; 793 } 794 if ((eyeval > -0.25) && (eyeval < 1.50) 795 && (diffeyecloudcat >= 2) && (eyefft <= 2) 796 && (cwtpart > 6.0) && (scenecloud <= 1) 797 && (cloudfft <= 4) && (lasttno12 >= 3.5)) { 798 sceneeye = 1; /* PINHOLE EYE CHECK */ 799 } 800 } else { 801 scenecloud = 4; /* SHEAR */ 802 } 803 } 804 } else { 805 scenecloud = 0; /* UNIFORM */ 806 if (embdd) { 807 scenecloud = 1; /* EMBEDDED CENTER */ 808 } 809 /* PINHOLE EYE TEST */ 810 /* 811 * printf("sceneeye=%d diffeyecloudcat=%d eyefft=%d cwtpart=%f scenecloud=%d cloudfft=%d lasttno12=%f\n" 812 * , 813 * sceneeye,diffeyecloudcat,eyefft,cwtpart,scenecloud,cloudfft,lasttno12 814 * ); 815 */ 816 if ((rmwsizeman > 0) && (rmwsizeman < 5.0)) { 817 sceneeye = 1; 818 } 819 if ((eyeval > -0.25) && (eyeval < 1.50) && (diffeyecloudcat >= 2) 820 && (eyefft <= 2) && (cwtpart > 6.0) && (scenecloud <= 1) 821 && (cloudfft <= 4) && (lasttno12 >= 3.5)) { 822 sceneeye = 1; /* PINHOLE EYE CHECK */ 823 } 824 } 825 if ((scenecloud <= 2) && (sceneeye == 3)) { 826 /* for CDO TESTS */ 827 for (ixx = 2; ixx <= 6; ixx++) { /* DG,MG,LG,B,W */ 828 tempval = (float) StormAODTInfo.ebd_v72[ixx] + 273.16f; 829 /* printf("xlat=%f xlon=%f tempval=%f\n",xlat,xlon,tempval); */ 830 float[] out = aodtv72_cdoshearcalc(xlat, xlon, tempval, 1, 831 areadata_v72); 832 /* 833 * printf("CDO : ixx=%d cdosize=%f cdosize/111=%f \n",ixx,cdosize 834 * ,cdosize/111.0); 835 */ 836 if (ixx == 2) { 837 eyecdosize = out[0]; 838 } 839 } 840 } 841 842 /* 843 * odtcurrent_v72IR.eyescene=sceneeye; 844 * odtcurrent_v72IR.cloudscene=scenecloud; 845 * odtcurrent_v72IR.eyesceneold=-1; odtcurrent_v72IR.cloudsceneold=-1; 846 * odtcurrent_v72IR.eyecdosize=eyecdosize; 847 * odtcurrent_v72IR.ringcb=cbring; odtcurrent_v72IR.ringcbval=cbringval; 848 * odtcurrent_v72IR.ringcbvalmax=cbringvalmax; 849 * odtcurrent_v72IR.ringcblatmax=cbringlatmax; 850 * odtcurrent_v72IR.ringcblonmax=cbringlonmax; 851 */ 852 float[] odt = { sceneeye, scenecloud, eyecdosize, cbring, cbringval, 853 cbringvalmax, cbringlatmax, cbringlonmax, odtcurrent_v72IRRMW }; 854 855 return odt; 856 857 } 858 859 /** 860 * Determine radius of maximum wind based upon Jim Kossin's regression based 861 * scheme. 862 * 863 * @param odtcurrent 864 * @param areadata 865 * 866 * @return Array of floats that contains radius of maximum wind distance 867 * and eye size radius (km). 868 */ 869 static float[] aodtv72_rmw(StormAODTInfo.IRData odtcurrent, 870 StormAODTInfo.DataGrid areadata) 871 /* 872 * Determine radius of maximum wind based upon Jim Kossin's regression based 873 * scheme Inputs : ix0 - element location of center point iy0 - line 874 * location of center point Outputs : rmw - radius of maximum wind distance 875 * eyesize - eye size radius (km) -1 = error/eyewall not found 0 = radius 876 * found 877 */ 878 { 879 int ixmin, ixmax, iymin, iymax; 880 int ixc, iyc, i, ix, iy, idx1 = 0, idy1 = 0, idx2 = 0, idy2 = 0; 881 float dx1, dx2, dy1, dy2, dav, xlat, xlon, xclat, xclon, xangle; 882 float tcrit = 228.0f, warm = 223.0f; 883 float rmw; 884 float eyesize; 885 886 /* 887 * calculate cursorx/cursory from numx/numy... values should be 888 * 0.5*numx/y 889 */ 890 ixc = areadata.numx / 2; 891 iyc = areadata.numy / 2; 892 ixmax = Math.min(areadata.numx, ixc + 320); 893 ixmin = Math.max(0, ixc - 320); 894 iymax = Math.min(areadata.numy, iyc + 240); 895 iymin = Math.max(0, iyc - 240); 896 897 if (odtcurrent.cloudt >= warm) { 898 tcrit = (odtcurrent.eyet + (2.0f * odtcurrent.cloudt)) / 3.0f; 899 } 900 901 rmw = -99.9f; 902 eyesize = -99.9f; 903 /* iterate five times */ 904 for (i = 0; i < 5; i++) { 905 ix = ixc; 906 while (areadata.temp[iyc][ix] > tcrit) { 907 ix = ix - 1; 908 if (ix == ixmin) { 909 return new float[] { -99.9f, -99.9f }; 910 } 911 } 912 idx1 = ix; 913 ix = ixc; 914 while (areadata.temp[iyc][ix] > tcrit) { 915 ix = ix + 1; 916 if (ix == ixmax) { 917 return new float[] { -99.9f, -99.9f }; 918 } 919 } 920 idx2 = ix; 921 iy = iyc; 922 while (areadata.temp[iy][ixc] > tcrit) { 923 iy = iy - 1; 924 if (iy == iymin) { 925 return new float[] { -99.9f, -99.9f }; 926 } 927 } 928 idy1 = iy; 929 iy = iyc; 930 while (areadata.temp[iy][ixc] > tcrit) { 931 iy = iy + 1; 932 if (iy == iymax) { 933 return new float[] { -99.9f, -99.9f }; 934 } 935 } 936 idy2 = iy; 937 ixc = (int) ((((float) (idx1 + idx2)) / 2.0)); 938 iyc = (int) ((((float) (idy1 + idy2)) / 2.0)); 939 } 940 xclat = areadata.lat[iyc][ixc]; 941 xclon = areadata.lon[iyc][ixc]; 942 xlat = areadata.lat[iyc][idx1]; 943 xlon = areadata.lon[iyc][idx1]; 944 float[] out = aodtv72_distance(xlat, xlon, xclat, xclon, 1); 945 dx1 = out[0]; 946 947 xlat = areadata.lat[iyc][idx2]; 948 xlon = areadata.lon[iyc][idx2]; 949 out = aodtv72_distance(xlat, xlon, xclat, xclon, 1); 950 dx2 = out[0]; 951 xlat = areadata.lat[idy1][ixc]; 952 xlon = areadata.lon[idy1][ixc]; 953 out = aodtv72_distance(xlat, xlon, xclat, xclon, 1); 954 dy1 = out[0]; 955 956 xlat = areadata.lat[idy2][ixc]; 957 xlon = areadata.lon[idy2][ixc]; 958 out = aodtv72_distance(xlat, xlon, xclat, xclon, 1); 959 dy2 = out[0]; 960 961 dav = (dx1 + dx2 + dy1 + dy2) / 4.0f; 962 if (dav > 0.0) { 963 rmw = 2.8068f + (0.8361f * dav); /* Howard's coeffs */ 964 eyesize = dav; 965 } 966 float[] outf = { rmw, eyesize }; 967 968 return outf; 969 } 970 971 /** 972 * Determine storm location using 10^ Log-spiral analysis. Algorithm will 973 * attempt to match the spiral with the image pixels at or below the 974 * threshold temperature based on BD-enhancement curve values 975 * 976 * @param inlat center latitude of analysis grid 977 * @param inlon center longitude of analysis grid 978 * @param searchtemp temperature threshold value 979 * @param searchtype 1=search at single point; 2=search over 2^box 980 * @param odtcurrent 981 * @param areadata 982 * 983 * @return Array of floats containing: {@code bestspiral} - number of 984 * consecutive arcs through which spiral passes, {@code bestlat} - best 985 * latitude location from analysis, and {@code bestlon} - best longitude 986 * location from analysis. 987 */ 988 static float[] aodtv72_logspiral(float inlat, float inlon, 989 float searchtemp, int searchtype, StormAODTInfo.IRData odtcurrent, 990 StormAODTInfo.DataGrid areadata) 991 /* 992 * Determine storm location using 10^ Log-spiral analysis. Algorithm will 993 * attempt to match the spiral with the image pixels at or below the 994 * threshold temperature based on BD-enhancement curve values Inputs : inlat 995 * - center latitude of analysis grid inlon - center longitude of analysis 996 * grid searchtemp - temperature threshold value searchtype - 1=search at 997 * single point;2=search over 2^box Outputs : bestlat - best latitude 998 * location from analysis bestlon - best longitude location from analysis 999 * bestspiral - number of consecutive arcs through which spiral passes 1000 */ 1001 { 1002 1003 int ixx, iyy, izz, iskip, rotfac = 0, theta, b; 1004 int imaxx, iminx, imaxy, iminy, ycount; 1005 int spiralconsec, maxconsec, spiralbest, spiralbestrot, bestrot = 0; 1006 float xrad = 57.29578f, A = 25.0f, B = 10.0f / xrad; 1007 float maxx, minx, maxy, miny, xmindist; 1008 float xlat = 0.f, xlon = 0.f; 1009 float ylatdiff, ylondiff, xres = 6.0f; 1010 float thetax, r, thetaval, thetaval2; 1011 /* float zlat[bufsiz],zlon[bufsiz]; */ 1012 float[] zlat, zlon, lres; 1013 int np, ipt; 1014 int bestspiral; 1015 float bestlat = 0.f; 1016 float bestlon = 0.f; 1017 1018 if (odtcurrent.sattype == 5) { 1019 xres = 12.0f; 1020 } 1021 /* allocate memory */ 1022 int nn = areadata.numx * areadata.numy; 1023 zlat = new float[nn]; 1024 zlon = new float[nn]; 1025 1026 if (searchtype == 2) { 1027 /* search over 2.0 degree box */ 1028 maxx = inlat + 1.0f; 1029 minx = inlat - 1.0f; 1030 maxy = inlon + 1.0f; 1031 miny = inlon - 1.0f; 1032 imaxx = (int) (maxx * 100.0); 1033 iminx = (int) (minx * 100.0); 1034 imaxy = (int) (maxy * 100.0); 1035 iminy = (int) (miny * 100.0); 1036 } else { 1037 /* search at a single point */ 1038 imaxx = (int) (inlat * 100.0); 1039 iminx = (int) (inlat * 100.0); 1040 imaxy = (int) (inlon * 100.0); 1041 iminy = (int) (inlon * 100.0); 1042 } 1043 1044 bestspiral = 0; 1045 1046 /* initialize arrays */ 1047 np = 0; 1048 for (ixx = 0; ixx < areadata.numx; ixx++) { 1049 for (iyy = 0; iyy < areadata.numy; iyy++) { 1050 if (areadata.temp[iyy][ixx] <= searchtemp) { 1051 zlat[np] = areadata.lat[iyy][ixx]; 1052 zlon[np] = areadata.lon[iyy][ixx]; 1053 np++; 1054 } 1055 } 1056 } 1057 1058 /* loop through x-axis/elements of analysis grid box */ 1059 for (ixx = iminx; ixx <= imaxx; ixx = ixx + 20) { 1060 xlat = (float) ixx / 100.0f; 1061 /* loop through y-axis/lines of analysis grid box */ 1062 for (iyy = iminy; iyy <= imaxy; iyy = iyy + 20) { 1063 xlon = (float) iyy / 100.0f; 1064 iskip = 0; 1065 /* determine distance from each point in box to current location */ 1066 if (searchtype == 2) { 1067 xmindist = 12.0f; 1068 for (izz = 0; izz < np; izz++) { 1069 float[] out3 = aodtv72_distance(xlat, xlon, zlat[izz], 1070 zlon[izz], 1); 1071 if (out3[0] <= xmindist) { 1072 /* 1073 * if the lat/lon point is too close to cold cloud 1074 * tops, do not calculate log spiral at this point. 1075 * Trying to eliminate "false" arc locations by 1076 * forcing the system to use some of the arc points 1077 * on the spiral away from the start of the spiral 1078 * (were getting "false echos" without this". 1079 */ 1080 iskip = 1; 1081 break; 1082 } 1083 } 1084 } 1085 1086 spiralbest = 0; 1087 spiralbestrot = 0; 1088 /* 1089 * if arc location passes analysis above, proceed with placement 1090 * of spiral 1091 */ 1092 if (iskip == 0) { 1093 /* 1094 * rotate the arc spiral thru entire revolution at 30 degree 1095 * intervals 1096 */ 1097 for (rotfac = 0; rotfac <= 330; rotfac = rotfac + 30) { 1098 spiralconsec = 0; 1099 maxconsec = 0; 1100 1101 /* 1102 * calculate position of each point on spiral from 0 to 1103 * 540^ 1104 */ 1105 for (theta = 0; theta <= 540; theta = theta + 15) { 1106 thetax = (float) theta / xrad; 1107 r = A * (float) Math.exp((B * thetax)); 1108 thetaval = (float) theta + (float) rotfac; 1109 if (xlat < 0.0) { 1110 thetaval = (float) (-1 * theta) 1111 + (float) rotfac; 1112 } 1113 thetaval2 = thetaval + 180.0f; 1114 float[] out2 = aodtv72_distance2(xlat, xlon, r, 1115 thetaval2); 1116 ycount = 0; 1117 for (izz = 0; izz < np; izz++) { 1118 ylatdiff = Math.abs(out2[0] - zlat[izz]); 1119 ylondiff = Math.abs(out2[1] - zlon[izz]); 1120 /* 1121 * if a point is within 0.1^ latitude/longitude 1122 * determine distance 1123 */ 1124 if ((ylatdiff <= 0.1) && (ylondiff <= 0.1)) { 1125 float[] out3 = aodtv72_distance(out2[0], 1126 out2[1], zlat[izz], zlon[izz], 1); 1127 /* 1128 * if distance from spiral point is within 1129 * 6km from an accepted temperature 1130 * threshold point, count it 1131 */ 1132 if (out3[0] <= xres) { 1133 ycount++; 1134 } 1135 } 1136 } 1137 /* 1138 * if there are 4 or more threshold points 1139 * associated with each spiral point, count within 1140 * consecutive spiral point counter 1141 */ 1142 if (ycount >= 4) { 1143 spiralconsec++; 1144 /* 1145 * save spiral that has the maximum consecutive 1146 * spiral counts for each rotation though 360^ 1147 * at each center location 1148 */ 1149 if (spiralconsec > maxconsec) { 1150 maxconsec = spiralconsec; 1151 } 1152 } else { 1153 spiralconsec = 0; 1154 } 1155 /* 1156 * if this spiral has the greatest number of 1157 * consecutive spiral points, save the location and 1158 * number of points 1159 */ 1160 if (maxconsec > spiralbest) { 1161 spiralbest = maxconsec; 1162 spiralbestrot = rotfac; 1163 } 1164 } 1165 } /* rotfac loop */ 1166 if (spiralbest > bestspiral) { 1167 bestspiral = spiralbest; 1168 bestlat = xlat; 1169 bestlon = xlon; 1170 bestrot = spiralbestrot; 1171 } 1172 } /* iskip if */ 1173 } /* iyy loop */ 1174 } /* ixx loop */ 1175 1176 /* free memory */ 1177 1178 /* load array for best spiral band */ 1179 ipt = 0; 1180 for (theta = 0; theta <= 540; theta = theta + 15) { 1181 thetax = (float) theta / xrad; 1182 r = A * (float) Math.exp((B * thetax)); 1183 thetaval = (float) theta + (float) bestrot; 1184 if (xlat < 0.0) { 1185 thetaval = (float) (-1 * theta) + (float) rotfac; 1186 } 1187 thetaval2 = thetaval + 180.0f; 1188 float[] out2 = aodtv72_distance2(bestlat, bestlon, r, thetaval2); 1189 /* load array for external plotting of spiral band */ 1190 // spiralband_v72[0][ipt]=out2[0]; 1191 // spiralband_v72[1][ipt]=out2[1]; 1192 ipt++; 1193 } 1194 1195 float[] out = { bestspiral, bestlat, bestlon }; 1196 return out; 1197 1198 } 1199 1200 /** 1201 * _more_ 1202 * 1203 * @param bin 1204 * _more_ 1205 * @param nbin 1206 * _more_ 1207 * 1208 * @return _more_ 1209 */ 1210 static float[] aodtv72_calcskew(float[] bin, int nbin) 1211 /* 1212 * Calculate average, standard deviation, and skew for a given data set. 1213 * Inputs : bin - data array nbin - number of points in data array Outputs : 1214 * ave - average value in array stdv - standard deviation of data array skew 1215 * - skew of data array histogram 1216 */ 1217 { 1218 int ixx; 1219 float xave, xstdv, xskew; 1220 float a2sum = 0.0f, a3sum = 0.0f, nsum = 0.0f; 1221 float ave, stdv, skew; 1222 1223 for (ixx = 0; ixx < nbin; ixx++) { 1224 nsum = nsum + bin[ixx]; 1225 } 1226 /* compute average value of data array */ 1227 xave = nsum / (float) nbin; 1228 1229 for (ixx = 0; ixx < nbin; ixx++) { 1230 a2sum = a2sum + ((bin[ixx] - xave) * (bin[ixx] - xave)); 1231 a3sum = a3sum 1232 + ((bin[ixx] - xave) * (bin[ixx] - xave) * (bin[ixx] - xave)); 1233 } 1234 if (nbin <= 1) { 1235 xstdv = 0.0f; 1236 xskew = 0.0f; 1237 } else { 1238 /* calculate standard deviation of data array */ 1239 xstdv = (float) Math.sqrt((1.0f / ((float) nbin - 1.0f)) * a2sum); 1240 /* calculate skew of data array */ 1241 xskew = (float) ((1.0f / ((float) nbin - 1.0f)) * a3sum) 1242 / (xstdv * xstdv * xstdv); 1243 } 1244 1245 /* return desired values */ 1246 ave = xave; 1247 stdv = xstdv; 1248 skew = xskew; 1249 float[] out = { ave, stdv, skew }; 1250 return out; 1251 } 1252 1253 /** 1254 * _more_ 1255 * 1256 * @param x1 1257 * _more_ 1258 * @param x2 1259 * _more_ 1260 * @param i2 1261 * _more_ 1262 * @param ndim 1263 * _more_ 1264 */ 1265 static void aodtv72_xxsort(float[] x1, float[] x2, float[] i2, int ndim) { 1266 int ih, i; 1267 float x, top; 1268 1269 ih = aodtv72_fminx(x1); 1270 x = x1[ih]; 1271 top = x - 1.0f; 1272 for (i = 0; i < ndim; i++) { 1273 ih = aodtv72_fmaxx(x1); 1274 i2[i] = ih; 1275 x2[i] = x1[ih]; 1276 x1[ih] = top; 1277 } 1278 } 1279 1280 /** 1281 * _more_ 1282 * 1283 * @param f 1284 * _more_ 1285 * 1286 * @return _more_ 1287 */ 1288 static int aodtv72_fminx(float[] f) { 1289 int i; 1290 int im; 1291 int ndim = f.length; 1292 float x; 1293 1294 im = 0; 1295 x = f[0]; 1296 for (i = 1; i < ndim; i++) { 1297 if (f[i] < x) { 1298 x = f[i]; 1299 im = i; 1300 } 1301 } 1302 return im; 1303 } 1304 1305 /** 1306 * _more_ 1307 * 1308 * @param f 1309 * _more_ 1310 * 1311 * @return _more_ 1312 */ 1313 static int aodtv72_fmaxx(float[] f) { 1314 int i; 1315 int im; 1316 int ndim = f.length; 1317 float x; 1318 1319 im = 0; 1320 x = f[0]; 1321 for (i = 1; i < ndim; i++) { 1322 if (f[i] > x) { 1323 x = f[i]; 1324 im = i; 1325 } 1326 } 1327 return im; 1328 } 1329 1330 /** 1331 * _more_ 1332 * 1333 * @param cbd 1334 * _more_ 1335 * 1336 * @return _more_ 1337 */ 1338 static float[] aodtv72_fft(float[] cbd) { 1339 int ixx, idxc, iok; 1340 double[] xr, xi, magn; 1341 double dx, a, x; 1342 float dx2; 1343 int idxcnt; 1344 1345 xr = new double[64]; 1346 xi = new double[64]; 1347 magn = new double[64]; 1348 /* initialize arrays */ 1349 for (ixx = 1; ixx <= 64; ++ixx) { 1350 xr[ixx - 1] = cbd[ixx - 1]; 1351 xi[ixx - 1] = (float) 0.; 1352 magn[ixx - 1] = (float) 0.; 1353 } 1354 /* call FFT routine */ 1355 iok = aodtv72_dfft(xr, xi, 64); 1356 if (iok <= 0) { 1357 return null; 1358 } 1359 for (ixx = 1; ixx <= 64; ++ixx) { 1360 magn[ixx - 1] = aodtv72_cmplx_abs(xr[ixx - 1], xi[ixx - 1]); 1361 } 1362 /* compute number of harmonics and magnitude of main harmonic */ 1363 dx = (float) 0.; 1364 idxc = 0; 1365 for (ixx = 2; ixx <= 31; ++ixx) { 1366 a = magn[ixx - 2]; 1367 x = magn[ixx - 1]; 1368 dx += (x + a) / (float) 2.; 1369 if ((magn[ixx - 1] > magn[ixx - 2]) && (magn[ixx - 1] > magn[ixx])) { 1370 ++idxc; 1371 } 1372 } 1373 dx2 = (float) (dx / magn[0]); 1374 1375 idxcnt = idxc; 1376 1377 float[] out = { dx2, idxcnt }; 1378 1379 return out; 1380 1381 } /* MAIN__ */ 1382 1383 /** 1384 * _more_ 1385 * 1386 * @param real 1387 * _more_ 1388 * @param imag 1389 * _more_ 1390 * 1391 * @return _more_ 1392 */ 1393 static double aodtv72_cmplx_abs(double real, double imag) { 1394 double temp; 1395 if (real < 0) { 1396 real = -real; 1397 } 1398 if (imag < 0) { 1399 imag = -imag; 1400 } 1401 if (imag > real) { 1402 temp = real; 1403 real = imag; 1404 imag = temp; 1405 } 1406 if ((real + imag) == real) { 1407 return (real); 1408 } 1409 1410 temp = imag / real; 1411 temp = real * Math.sqrt(1.0 + temp * temp); /* overflow!! */ 1412 return (temp); 1413 1414 } 1415 1416 /* 1417 * A Duhamel-Hollman split-radix dif fft Ref: Electronics Letters, Jan. 5, 1418 * 1984 Complex input and output data in arrays x and y Length is n. 1419 */ 1420 1421 /** 1422 * _more_ 1423 * 1424 * @param x 1425 * _more_ 1426 * @param y 1427 * _more_ 1428 * @param np 1429 * _more_ 1430 * 1431 * @return _more_ 1432 */ 1433 static int aodtv72_dfft(double[] x, double[] y, int np) { 1434 1435 int i, j, k, m, n, i0, i1, i2, i3, is, id, n1, n2, n4; 1436 double a, e, a3, cc1, ss1, cc3, ss3, r1, r2, s1, s2, s3, xt; 1437 1438 double px[] = new double[x.length + 1]; 1439 double py[] = new double[y.length + 1]; 1440 1441 px[0] = 0; 1442 py[0] = 0; 1443 for (i = 0; i < x.length; i++) { 1444 px[i + 1] = x[i]; 1445 py[i + 1] = y[i]; 1446 } 1447 1448 i = 2; 1449 m = 1; 1450 1451 while (i < np) { 1452 i = i + i; 1453 m = m + 1; 1454 } 1455 ; 1456 n = i; 1457 if (n != np) { 1458 for (i = np + 1; i <= n; i++) { 1459 px[i] = 0.0; 1460 py[i] = 0.0; 1461 } 1462 /* printf("\nuse %d point fft",n); */ 1463 } 1464 1465 n2 = n + n; 1466 for (k = 1; k <= m - 1; k++) { 1467 n2 = n2 / 2; 1468 n4 = n2 / 4; 1469 e = 2.0 * PI / n2; 1470 a = 0.0; 1471 for (j = 1; j <= n4; j++) { 1472 a3 = 3.0 * a; 1473 cc1 = Math.cos(a); 1474 ss1 = Math.sin(a); 1475 cc3 = Math.cos(a3); 1476 ss3 = Math.sin(a3); 1477 a = j * e; 1478 is = j; 1479 id = 2 * n2; 1480 while (is < n) { 1481 for (i0 = is; i0 <= n - 1; i0 = i0 + id) { 1482 i1 = i0 + n4; 1483 i2 = i1 + n4; 1484 i3 = i2 + n4; 1485 r1 = px[i0] - px[i2]; 1486 px[i0] = px[i0] + px[i2]; 1487 r2 = px[i1] - px[i3]; 1488 px[i1] = px[i1] + px[i3]; 1489 s1 = py[i0] - py[i2]; 1490 py[i0] = py[i0] + py[i2]; 1491 s2 = py[i1] - py[i3]; 1492 py[i1] = py[i1] + py[i3]; 1493 s3 = r1 - s2; 1494 r1 = r1 + s2; 1495 s2 = r2 - s1; 1496 r2 = r2 + s1; 1497 px[i2] = r1 * cc1 - s2 * ss1; 1498 py[i2] = -s2 * cc1 - r1 * ss1; 1499 px[i3] = s3 * cc3 + r2 * ss3; 1500 py[i3] = r2 * cc3 - s3 * ss3; 1501 } 1502 is = 2 * id - n2 + j; 1503 id = 4 * id; 1504 } 1505 } 1506 } 1507 1508 /* 1509 * ---------------------Last stage, length=2 1510 * butterfly--------------------- 1511 */ 1512 is = 1; 1513 id = 4; 1514 while (is < n) { 1515 for (i0 = is; i0 <= n; i0 = i0 + id) { 1516 i1 = i0 + 1; 1517 r1 = px[i0]; 1518 px[i0] = r1 + px[i1]; 1519 px[i1] = r1 - px[i1]; 1520 r1 = py[i0]; 1521 py[i0] = r1 + py[i1]; 1522 py[i1] = r1 - py[i1]; 1523 } 1524 is = 2 * id - 1; 1525 id = 4 * id; 1526 } 1527 1528 /* 1529 * c--------------------------Bit reverse counter 1530 */ 1531 j = 1; 1532 n1 = n - 1; 1533 for (i = 1; i <= n1; i++) { 1534 if (i < j) { 1535 xt = px[j]; 1536 px[j] = px[i]; 1537 px[i] = xt; 1538 xt = py[j]; 1539 py[j] = py[i]; 1540 py[i] = xt; 1541 } 1542 k = n / 2; 1543 while (k < j) { 1544 j = j - k; 1545 k = k / 2; 1546 } 1547 j = j + k; 1548 } 1549 1550 /* 1551 * for (i = 1; i<=16; i++) printf("%d %g %gn",i,*(px+i),(py+i)); 1552 */ 1553 for (i = 0; i < x.length - 1; i++) { 1554 x[i] = px[i + 1]; 1555 y[i] = py[i + 1]; 1556 } 1557 return (n); 1558 } 1559 1560 /** 1561 * _more_ 1562 * 1563 * @param keyer 1564 * _more_ 1565 * @param areadata 1566 * _more_ 1567 * 1568 * @return _more_ 1569 */ 1570 static public StormAODTInfo.IRData aodtv72_gettemps(int keyer, 1571 StormAODTInfo.DataGrid areadata) 1572 /* 1573 * Determine eye and coldest-warmest cloud temperatures. Inputs : keyer - 1574 * eye radius (km) Outputs : odtcurrent_v72 - structure containing current 1575 * image information Return : -51 : eye, cloud, or warmest temperature 1576 * <-100C or >+40C 1577 */ 1578 { 1579 1580 int cursorx, cursory; 1581 int iok, idist, iret; 1582 float cursortemp, eyet, cloudt; 1583 float cursorlat, cursorlon, warmlat, warmlon, eyeangle, eyedist; 1584 StormAODTInfo sinfo = new StormAODTInfo(); 1585 StormAODTInfo.IRData odtcurrent = new StormAODTInfo.IRData(); 1586 1587 /* 1588 * calculate cursorx/cursory from numx/numy... values should be 1589 * 0.5*numx/y 1590 */ 1591 cursorx = areadata.numx / 2; 1592 cursory = areadata.numy / 2; 1593 cursortemp = areadata.temp[cursory][cursorx]; 1594 cursorlat = areadata.lat[cursory][cursorx]; 1595 cursorlon = areadata.lon[cursory][cursorx]; 1596 1597 // tcircfirst_v72=(struct ringdata *)calloc(1,sizeof(struct ringdata)); 1598 1599 /* load array containing data on rings around center location */ 1600 List<StormAODTInfo.RingData> tcircfirst = aodtv72_readcirc(cursorx, 1601 cursory, areadata); 1602 1603 iret = 0; 1604 /* compute eye/warmest pixel temperature */ 1605 StormAODTInfo.RingData eye = aodtv72_calceyetemp(keyer, cursortemp, 1606 tcircfirst); 1607 eyet = eye.temp; 1608 if ((eyet < -100.0) || (eyet > 40.0)) { 1609 iret = -51; 1610 } 1611 1612 if (keyer == StormAODTInfo.keyerA_v72) { 1613 /* set warmest pixel temperature */ 1614 odtcurrent.warmt = eyet; 1615 /* calculate warmest pixel location */ 1616 float[] out2 = aodtv72_distance2(cursorlat, cursorlon, eye.dist, 1617 eye.angle); 1618 1619 odtcurrent.warmlatitude = out2[0]; // warmlat; 1620 odtcurrent.warmlongitude = out2[1]; // warmlon; 1621 /* store forecast location temperature in eyet */ 1622 odtcurrent.eyet = cursortemp - 273.16f; 1623 // free(tcircfirst_v72); 1624 } else { 1625 /* set eye temperature */ 1626 odtcurrent.eyet = eyet; 1627 /* set cloud temperature and ring distance */ 1628 float[] out3 = aodtv72_calccloudtemp(tcircfirst); 1629 cloudt = out3[0]; 1630 idist = (int) out3[1]; 1631 if ((eyet < -100.0) || (eyet > 40.0)) { 1632 iret = -51; 1633 } 1634 odtcurrent.cwcloudt = cloudt; 1635 odtcurrent.cwring = idist; 1636 } 1637 if (iret == -51) 1638 return null; 1639 return odtcurrent; 1640 } 1641 1642 /** 1643 * _more_ 1644 * 1645 * @param keyer 1646 * _more_ 1647 * @param cursortemp 1648 * _more_ 1649 * @param tcircfirst 1650 * _more_ 1651 * 1652 * @return _more_ 1653 */ 1654 static public StormAODTInfo.RingData aodtv72_calceyetemp(int keyer, 1655 float cursortemp, List<StormAODTInfo.RingData> tcircfirst) 1656 /* 1657 * Determine eye/warmest temperature by examining the satellite data between 1658 * 0 and 24/75 km from the storm center location. Eye temperature will be 1659 * warmest temperature found. Inputs : keyer - analysis region radius 1660 * cursortemp - temperature of pixel at cursor location Outputs : eyeangle - 1661 * angle to warmest temperature in eye region (if found) eyedist - distance 1662 * to warmest temperature in eye region (if found) Return : return value is 1663 * warmest eye temperature 1664 */ 1665 { 1666 1667 List<StormAODTInfo.RingData> tcirc; 1668 StormAODTInfo.RingData eye = new StormAODTInfo.RingData(); 1669 /* set eye temp to cursor location temperature */ 1670 eye.temp = cursortemp; 1671 1672 eye.dist = 0.0f; 1673 eye.angle = 0.0f; 1674 tcirc = tcircfirst; 1675 Iterator it = tcirc.iterator(); 1676 while (it.hasNext()) { 1677 StormAODTInfo.RingData rd = (StormAODTInfo.RingData) it.next(); 1678 if (rd.dist <= (float) keyer) { 1679 if (rd.temp > eye.temp) { 1680 eye.temp = rd.temp; 1681 eye.dist = rd.dist; 1682 eye.angle = rd.angle; 1683 } 1684 } 1685 1686 } 1687 1688 /* convert temperature to C from K */ 1689 eye.temp = (eye.temp - 273.16f); 1690 1691 return eye; 1692 } 1693 1694 /** 1695 * _more_ 1696 * 1697 * @param tcircfirst 1698 * _more_ 1699 * 1700 * @return _more_ 1701 */ 1702 static public float[] aodtv72_calccloudtemp( 1703 List<StormAODTInfo.RingData> tcircfirst) 1704 /* 1705 * Determine surrounding cloud top region temperature by examining the 1706 * satellite data between kstart_v72(24 km) and kend_v72 (136 km) from the 1707 * storm center location. Cloud temperature will be the coldest value of an 1708 * array of warmest ring temperatures (4 km ring sizes for 4 km infrared 1709 * data). This is the "coldest-warmest" pixel. Inputs : none Outputs : 1710 * irdist - distance (km) from center to cloud top temperature value 1711 * (distance to ring) return value is cloud top temperature value 1712 */ 1713 { 1714 int iyy, b; 1715 int numrings = (StormAODTInfo.kend_v72 - StormAODTInfo.kstart_v72) 1716 / StormAODTInfo.kres_v72; 1717 int kring; 1718 /* float maxtemp[200]; */ 1719 float[] maxtemp; 1720 1721 List<StormAODTInfo.RingData> tcirc; 1722 1723 float cloudtemp = 10000.0f; 1724 1725 float[] out = new float[2]; 1726 int irdist = 0; 1727 /* initialize maxtemp array */ 1728 maxtemp = new float[numrings]; 1729 1730 tcirc = tcircfirst; 1731 Iterator it = tcirc.iterator(); 1732 while (it.hasNext()) { 1733 StormAODTInfo.RingData rd = (StormAODTInfo.RingData) it.next(); 1734 if ((rd.dist >= (float) StormAODTInfo.kstart_v72) 1735 && (rd.dist < (float) StormAODTInfo.kend_v72)) { 1736 kring = ((int) rd.dist - StormAODTInfo.kstart_v72) 1737 / StormAODTInfo.kres_v72; 1738 if (rd.temp > maxtemp[kring]) { 1739 maxtemp[kring] = rd.temp; 1740 } 1741 } 1742 } 1743 1744 /* search maxtemp array for coldest temperature */ 1745 for (iyy = 0; iyy < numrings; iyy++) { 1746 if ((maxtemp[iyy] < cloudtemp) && (maxtemp[iyy] > 160.0)) { 1747 cloudtemp = maxtemp[iyy]; 1748 irdist = (iyy * StormAODTInfo.kres_v72) 1749 + StormAODTInfo.kstart_v72; 1750 } 1751 } 1752 1753 cloudtemp = (cloudtemp - 273.16f); 1754 out[0] = cloudtemp; 1755 out[1] = irdist; 1756 return out; 1757 } 1758 1759 /** 1760 * _more_ 1761 * 1762 * @param ixc 1763 * _more_ 1764 * @param iyc 1765 * _more_ 1766 * @param areadata 1767 * _more_ 1768 * 1769 * @return _more_ 1770 */ 1771 static public List<StormAODTInfo.RingData> aodtv72_readcirc(int ixc, 1772 int iyc, StormAODTInfo.DataGrid areadata) 1773 /* 1774 * Read array of satellite data temperatures and load array containing 1775 * temperatures and lat/lon positions. Inputs : ixc - element location of 1776 * center point iyc - line location of center point Outputs : global 1777 * structure tcirc containing temperature/locations 1778 */ 1779 { 1780 int ixx, iyy, b, count = 0; 1781 float xclat, xclon, xlat, xlon, xdist, xangle; 1782 StormAODTInfo sinfo = new StormAODTInfo(); 1783 StormAODTInfo.RingData tcirc = null; 1784 1785 List<StormAODTInfo.RingData> tcircfirst = new ArrayList(); 1786 1787 /* obtain center/cursor x,y position */ 1788 xclat = areadata.lat[iyc][ixc]; 1789 xclon = areadata.lon[iyc][ixc]; 1790 /* load tcirc array with distance, angle, and temp info */ 1791 1792 for (iyy = 0; iyy < areadata.numy; iyy++) { 1793 for (ixx = 0; ixx < areadata.numx; ixx++) { 1794 xlat = areadata.lat[iyy][ixx]; 1795 xlon = areadata.lon[iyy][ixx]; 1796 float[] outDA = aodtv72_distance(xlat, xlon, xclat, xclon, 1); 1797 xdist = outDA[0]; 1798 xangle = outDA[1]; 1799 if (xdist <= (float) (StormAODTInfo.kend_v72 + 80)) { /* 1800 * add 80.0 1801 * to allow 1802 * for 1803 * correct 1804 * calculation 1805 * of 1806 * annulus 1807 * temp 1808 */ 1809 count++; 1810 tcirc = new StormAODTInfo.RingData(xdist, xangle, 1811 areadata.temp[iyy][ixx]); 1812 tcircfirst.add(tcirc); 1813 // tcirc.nextrec=NULL; /* make pointer for last record equal 1814 // to 0 */ 1815 } 1816 } 1817 } 1818 1819 return tcircfirst; 1820 } 1821 1822 /** 1823 * _more_ 1824 * 1825 * @param rrlat 1826 * _more_ 1827 * @param rrlon 1828 * _more_ 1829 * @param pplat 1830 * _more_ 1831 * @param pplon 1832 * _more_ 1833 * @param iunit 1834 * _more_ 1835 * 1836 * @return _more_ 1837 */ 1838 static float[] aodtv72_distance(float rrlat, float rrlon, float pplat, 1839 float pplon, int iunit) 1840 /* 1841 * Calculate distance and angle between two points (rrlat,rrlon and 1842 * pplat,pplon). Inputs : rrlat - latitude of starting point rrlon - 1843 * longitude of starting point rrlat - latitude of ending point rrlon - 1844 * longitude of ending point iunit - flag for output unit type 1845 * (1-km,2-mi,3-nmi) Outputs : dist - distance between two points ang - 1846 * angle between two points (0=north) 1847 */ 1848 { 1849 float z = 0.017453292f, r = 6371.0f; 1850 float rlat, rlon, plat, plon; 1851 float crlat, crlon, srlat, srlon; 1852 float cplat, cplon, splat, splon; 1853 float aa, xx, yy, zz, idist, xdist, xang; 1854 float[] out = new float[2]; 1855 float dist; 1856 float ang; 1857 1858 idist = 0.0f; 1859 rlat = rrlat * z; 1860 rlon = rrlon * z; 1861 plat = pplat * z; 1862 plon = pplon * z; 1863 crlat = (float) Math.cos(rlat); 1864 crlon = (float) Math.cos(rlon); 1865 srlat = (float) Math.sin(rlat); 1866 srlon = (float) Math.sin(rlon); 1867 cplat = (float) Math.cos(plat); 1868 cplon = (float) Math.cos(plon); 1869 splat = (float) Math.sin(plat); 1870 splon = (float) Math.sin(plon); 1871 xx = (cplat * cplon) - (crlat * crlon); 1872 yy = (cplat * splon) - (crlat * srlon); 1873 zz = splat - srlat; 1874 aa = (xx * xx) + (yy * yy) + (zz * zz); 1875 idist = (float) Math.sqrt(aa); 1876 /* xdist is distance in kilometers */ 1877 xdist = 2.0f * (float) Math.asin(idist / 2.0f) * r; 1878 1879 if (iunit == 2) { 1880 xdist = ((69.0f * xdist) + 55f) / 111.0f; /* conversion to miles */ 1881 } 1882 if (iunit == 3) { 1883 xdist = ((60.0f * xdist) + 55f) / 111.0f; /* 1884 * conversion to nautical 1885 * miles 1886 */ 1887 } 1888 dist = xdist; 1889 1890 xang = 0.0f; 1891 if (Math.abs(xdist) > 0.0001) { 1892 xang = (float) ((Math.sin(rlon - plon) * Math.sin((3.14159 / 2.0) 1893 - plat)) / Math.sin(idist)); 1894 } 1895 if (Math.abs(xang) > 1.0) { 1896 xang = (xang >= 0) ? 1 : -1; 1897 } 1898 xang = (float) Math.asin(xang) / z; 1899 if (plat < rlat) { 1900 xang = 180.0f - xang; 1901 } 1902 if (xang < 0.0) { 1903 xang = 360.0f + xang; 1904 } 1905 ang = xang; 1906 out[0] = dist; 1907 out[1] = ang; 1908 1909 return out; 1910 } 1911 1912 /** 1913 * _more_ 1914 * 1915 * @param rlat 1916 * _more_ 1917 * @param rlon 1918 * _more_ 1919 * @param xdist 1920 * _more_ 1921 * @param xang 1922 * _more_ 1923 * 1924 * @return _more_ 1925 */ 1926 static float[] aodtv72_distance2(float rlat, float rlon, float xdist, 1927 float xang) 1928 /* 1929 * Calculate a latitude and longitude position from an initial 1930 * latitude/longitude and distance/angle values. Inputs : rlat - initial 1931 * latitude rlon - initial longitude dist - distance from initial position 1932 * ang - angle from initial position Outputs : plat - derived latitude plon 1933 * - derived longitude 1934 */ 1935 { 1936 float z = 0.017453292f, z90 = 1.570797f; 1937 float clat, clon, cltv, cdis, cang; 1938 float qlat, qlon, argu, tv; 1939 int iang; 1940 float plat, plon; 1941 float[] out = new float[2]; 1942 1943 clat = (90.0f - rlat) * z; 1944 cltv = clat; 1945 clon = rlon * z; 1946 if (rlat < 0.0) { 1947 cltv = -(90.0f + rlat) * z; 1948 clon = (rlon - 180.0f) * z; 1949 xang = 360 - xang; 1950 } 1951 iang = (int) xang; 1952 cang = -(float) (((540 - iang) % 360) * z); 1953 cdis = (xdist / 111.1f) * z; 1954 qlat = (float) Math.acos((Math.cos(clat) * Math.cos(cdis)) 1955 + (Math.sin(clat) * Math.sin(cdis) * Math.cos(cang))); 1956 if (Math.abs(qlat) < 0.0000001) { 1957 qlon = 0.0f; 1958 } else { 1959 argu = (float) (Math.sin(cdis) * Math.sin(cang)) 1960 / (float) Math.sin(qlat); 1961 if (Math.abs(argu) > 1.0) { 1962 argu = (argu >= 0) ? 1 : -1; 1963 ; 1964 } 1965 qlon = (float) Math.asin(argu); 1966 tv = (float) Math.atan(Math.sin(z90 - cang)) 1967 / (float) Math.tan(z90 - cdis); 1968 if (tv > cltv) { 1969 qlon = (2.0f * z90) - qlon; 1970 } 1971 } 1972 qlon = clon - qlon; 1973 plat = 90.0f - (qlat / z); 1974 plon = (float) ((int) (10000 * (qlon / z)) % 3600000) / 10000.0f; 1975 if (plon < -180) { 1976 plon = plon + 360; 1977 } 1978 1979 out[0] = plat; 1980 out[1] = plon; 1981 1982 return out; 1983 } 1984 1985 /** 1986 * _more_ 1987 * 1988 * @param date 1989 * _more_ 1990 * @param time 1991 * _more_ 1992 * 1993 * @return _more_ 1994 */ 1995 static double aodtv72_calctime(int date, int time) 1996 /* 1997 * Compute time in xxxxx.yyy units, where xxxxx is the day and yyy is the 1998 * percentage of the day. This routine will also correct for Y2K problems. 1999 * Inputs : date - Julian date time - time in HHMMSS format Outputs : 2000 * function return value 2001 */ 2002 { 2003 int iyy; 2004 float sec, min, hour, partday; 2005 double timeout; 2006 2007 if ((date % 1000) == 0) { 2008 return 0; 2009 } 2010 if (time < 0) { 2011 return 0; 2012 } 2013 2014 iyy = date / 1000; /* obtain year */ 2015 /* 2016 * check for millenium designation in the year. if it is not there, add 2017 * it onto the beginning 2018 */ 2019 if (iyy < 1900) { 2020 if (iyy > 70) { 2021 date = 1900000 + date; 2022 } else { 2023 date = 2000000 + date; 2024 } 2025 } 2026 2027 sec = ((float) (time % 100)) / 3600.0f; 2028 min = ((float) ((time / 100) % 100)) / 60.0f; 2029 hour = (float) (time / 10000); 2030 partday = (hour + min + sec) / 24.0f; 2031 timeout = (double) date + (double) partday; 2032 2033 return timeout; 2034 } 2035 2036 /** 2037 * Determine eye size or shear distance for a given scene. 2038 * 2039 * @param xlat Center latitude of analysis grid. 2040 * @param xlon Center longitude of analysis grid. 2041 * @param tempval Temperature threshold value to be used. 2042 * @param atype Analysis type (1: cdo size; 2: eye size; 3: shear distance). 2043 * @param areadata _more_ 2044 * 2045 * @return Array of floats containing: {@code valb} - eye/cdo radius 2046 * or shear distance, and {@code valc} - eye/cdo symmetry value or 0. 2047 */ 2048 static float[] aodtv72_cdoshearcalc(float xlat, float xlon, float tempval, 2049 int atype, StormAODTInfo.DataGrid areadata) 2050 /* 2051 * Determine eye size or shear distance for a given scene. Inputs : xlat - 2052 * center latitude of analysis grid xlon - center longitude of analysis grid 2053 * tempval - temperature threshold value to be used atype - analysis type 2054 * (1-cdo size,2-eye size,3-shear distance) Outputs : valb - eye/cdo radius 2055 * or shear distance valc - eye/cdo symmetry value or 0 2056 */ 2057 { 2058 2059 int ixx, iyy, np, b, numvalid = 0; 2060 float smalldist, vc, maxdist; 2061 float[] zlat, zlon; 2062 float a1, a2, a3, a4, v1, v2, v3; 2063 float valb = 0.f, valc = 0.f; 2064 /* allocate memory */ 2065 2066 np = 0; 2067 int nn = areadata.numx * areadata.numy; 2068 zlat = new float[nn]; 2069 zlon = new float[nn]; 2070 /* CDO size determination - RETURNS RADIUS */ 2071 if (atype == 1) { 2072 /* a1=999.9;a2=999.9;a3=999.9;a4=999.9; */ 2073 a1 = 300.0f; 2074 a2 = 300.0f; 2075 a3 = 300.0f; 2076 a4 = 300.0f; 2077 for (ixx = 0; ixx < areadata.numx; ixx++) { 2078 for (iyy = 0; iyy < areadata.numy; iyy++) { 2079 if (areadata.temp[iyy][ixx] > tempval) { 2080 zlat[np] = areadata.lat[iyy][ixx]; 2081 zlon[np] = areadata.lon[iyy][ixx]; 2082 np++; 2083 } 2084 } 2085 } 2086 /* 2087 * printf("np=%d numx*numy=%d\n",np,areadata_v72->numx*areadata_v72- 2088 * >numy); 2089 */ 2090 maxdist = 0.0f; 2091 if (np < (areadata.numx * areadata.numy)) { 2092 for (ixx = 0; ixx < np; ixx++) { 2093 float[] out = aodtv72_distance(xlat, xlon, zlat[ixx], 2094 zlon[ixx], 1); 2095 if (out[0] > maxdist) { 2096 maxdist = out[0]; 2097 } 2098 /* determine size of CDO */ 2099 if (out[0] > StormAODTInfo.keyerM_v72) { 2100 if ((Math.abs(out[1] - 45.0) <= 15.0) && (out[0] < a1)) { 2101 a1 = out[0]; 2102 } 2103 if ((Math.abs(out[1] - 135.0) <= 15.0) && (out[0] < a2)) { 2104 a2 = out[0]; 2105 } 2106 if ((Math.abs(out[1] - 225.0) <= 15.0) && (out[0] < a3)) { 2107 a3 = out[0]; 2108 } 2109 if ((Math.abs(out[1] - 315.0) <= 15.0) && (out[0] < a4)) { 2110 a4 = out[0]; 2111 } 2112 } 2113 } 2114 /* printf("a1=%f a2=%f a3=%f a4=%f\n",a1,a2,a3,a4); */ 2115 numvalid = 4; 2116 v3 = StormAODTInfo.keyerM_v72 + StormAODTInfo.kres_v72; 2117 if (a1 < v3) { 2118 numvalid--; 2119 } 2120 if (a2 < v3) { 2121 numvalid--; 2122 } 2123 if (a3 < v3) { 2124 numvalid--; 2125 } 2126 if (a4 < v3) { 2127 numvalid--; 2128 } 2129 } else { 2130 a1 = 0.0f; 2131 a2 = 0.0f; 2132 a3 = 0.0f; 2133 a4 = 0.0f; 2134 } 2135 if (numvalid < 3) { 2136 valb = 0.0f; 2137 valc = 0.0f; 2138 } else { 2139 a1 = Math.min(a1, maxdist); 2140 a2 = Math.min(a2, maxdist); 2141 a3 = Math.min(a3, maxdist); 2142 a4 = Math.min(a4, maxdist); 2143 valb = (a1 + a2 + a3 + a4) / 4.0f; 2144 /* *valc=A_ABS(((a1+a3)/2.0)-((a2+a4)/2.0))/2.0; */ 2145 v1 = a1 + a3; 2146 v2 = a2 + a4; 2147 vc = v1 / v2; 2148 vc = Math.max(vc, 1.0f / vc); 2149 valc = vc; 2150 } 2151 /* 2152 * printf("\nnp=%5.5d a1=%5.1f a2=%5.1f a3=%5.1f a4=%5.1f ",np,a1,a2 2153 * ,a3,a4); 2154 */ 2155 } 2156 2157 /* shear distance determination */ 2158 if (atype == 3) { 2159 a1 = 999.9f; 2160 a2 = 999.9f; 2161 a3 = 999.9f; 2162 a4 = 999.9f; 2163 for (ixx = 0; ixx < areadata.numx; ixx++) { 2164 for (iyy = 0; iyy < areadata.numy; iyy++) { 2165 if (areadata.temp[iyy][ixx] <= tempval) { 2166 zlat[np] = areadata.lat[iyy][ixx]; 2167 zlon[np] = areadata.lon[iyy][ixx]; 2168 np++; 2169 } 2170 } 2171 } 2172 smalldist = 999.9f; 2173 for (ixx = 0; ixx < np; ixx++) { 2174 float[] out = aodtv72_distance(xlat, xlon, zlat[ixx], 2175 zlon[ixx], 1); 2176 if (out[0] < smalldist) { 2177 smalldist = out[0]; 2178 } 2179 } 2180 valb = smalldist; 2181 valc = 0.0f; 2182 2183 } 2184 2185 /* free memory */ 2186 2187 float[] out = { valb, valc }; 2188 return out; 2189 2190 } 2191 2192}