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