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.data.cyclone;
030    
031    import java.util.List;
032    
033    import ucar.unidata.data.grid.GridUtil;
034    import ucar.unidata.util.IOUtil;
035    import ucar.unidata.util.StringUtil;
036    import visad.FlatField;
037    
038    /**
039     * Created by IntelliJ IDEA. User: yuanho Date: Feb 20, 2009 Time: 3:09:14 PM To
040     * change this template use File | Settings | File Templates.
041     */
042    
043    public class StormAODT {
044    
045            /** _more_ */
046            StormAODTInfo.IRData odtcurrent_v72IR;
047    
048            /** _more_ */
049            StormAODTInfo.DataGrid areadata_v72;
050    
051            /** _more_ */
052            boolean lauto = false;
053    
054            /** _more_ */
055            int idomain_v72, ixdomain_v72, ifixtype_v72, rmwsizeman_v72;
056    
057            /** _more_ */
058            int oland_v72;
059    
060            /** _more_ */
061            boolean osearch_v72;
062    
063            /** _more_ */
064            int ostartstr_v72;
065    
066            /** _more_ */
067            float osstr_v72;
068    
069            /** _more_ */
070            static double c1 = 1.191066E-5;
071    
072            /** _more_ */
073            static double c2 = 1.438833;
074    
075            public StormAODT() {
076    
077            }
078    
079            /**
080             * _more_
081             * 
082             * 
083             * @param satgrid
084             *            _more_
085             * @param cenlat
086             *            _more_
087             * @param cenlon
088             *            _more_
089             * @param posm
090             *            _more_
091             * @param curdate
092             *            _more_
093             * @param cursat
094             *            _more_
095             * @param g_domain
096             *            _more_
097             * 
098             * @return _more_
099             */
100            public StormAODTInfo.IRData aodtv72_drive(FlatField satgrid, float cenlat,
101                            float cenlon, int posm, double curdate, int cursat,
102                            String g_domain, int satId, int satChannel, boolean isTemperature) {
103    
104                    float ftmps, flats, flons, cenlon2;
105    
106                    int radius, irad, np, ii, jj, length;
107                    int idomain = 0;
108    
109                    /*
110                     * Set miscoptions flags in AODT
111                     */
112    
113                    int eyeSize = -99;
114                    oland_v72 = 0; /* allow AODT operation over land */
115                    osearch_v72 = false; /* search for maximum curved band position */
116                    rmwsizeman_v72 = eyeSize; /* eye size parameter */
117                    odtcurrent_v72IR = new StormAODTInfo.IRData();
118                    odtcurrent_v72IR.domain = idomain_v72;
119    
120                    /*
121                     * Set initial classification flag and value in AODT
122                     */
123    
124                    ostartstr_v72 = 0; /* user defined initial classification flag */
125                    osstr_v72 = 0.0f; /* starting initial classification value */
126    
127                    /*
128                     * Set image date/time info in AODT
129                     */
130    
131                    int iaodt = aodtv72_setIRimageinfo(curdate, cursat);
132    
133                    /*
134                     * Get storm center lat/lon
135                     */
136                    if (lauto == true) {
137                            // aodtv72_runautomode( nauto, fauto, imagefile, &cenlat, &cenlon,
138                            // &posm );
139                    }
140    
141                    /*
142                     * Set center location in AODT positioning method (1=interpolation,
143                     * 4=extrapolation, 0=error)
144                     */
145                    // posm = 1;
146                    aodtv72_setlocation(cenlat, cenlon, posm);
147    
148                    /*
149                     * Set domain FLAG in AODT
150                     */
151    
152                    if (g_domain.equalsIgnoreCase("AUTO")) {
153                            idomain = 0;
154                    }
155                    if (g_domain.equalsIgnoreCase("Atlantic")) {
156                            idomain = 1;
157                    }
158                    if (g_domain.equalsIgnoreCase("Pacific")) {
159                            idomain = 2;
160                    }
161                    if (g_domain.equalsIgnoreCase("Indian")) {
162                            idomain = 2;
163                    }
164    
165                    iaodt = aodtv72_setdomain(idomain);
166    
167                    /*
168                     * Retrieve temperatures from image. This to be done in IDV
169                     */
170    
171                    GridUtil.Grid2D g2d = null;
172                    float[][] temps = null;
173                    float[][][] satimage = null;
174                    float[][] lons = null;
175                    float[][] lats = null;
176                    int numx = 123;
177                    int numy = 123;
178    
179                    try {
180                            g2d = GridUtil.makeGrid2D(satgrid);
181                            lons = g2d.getlons();
182                            lats = g2d.getlats();
183    
184                    } catch (Exception re) {
185                    }
186    
187                    /* now spatial subset numx by numy */
188                    GridUtil.Grid2D g2d1 = spatialSubset(g2d, cenlat, cenlon, numx, numy);
189    
190                    satimage = g2d1.getvalues();
191                    float[][] temp0 = satimage[0];
192                    int imsorc = satId, imtype = satChannel;
193    
194                    if (isTemperature)
195                            temps = temp0;
196                    else
197                            temps = im_gvtota(numx, numy, temp0, imsorc, imtype);
198    
199                    /*
200                     * Load the IR image information in AODT init areadata_v72
201                     */
202    
203                    aodtv72_loadIRimage(temps, g2d1.getlats(), g2d1.getlons(), numx, numy);
204    
205                    /*
206                     * Set eye and cloud temperature values in AODT, return position for IR
207                     * image data read
208                     */
209    
210                    StormAODTInfo.IRData tvIR = aodtv72_seteyecloudtemp(
211                                    StormAODTInfo.keyerM_v72, areadata_v72);
212    
213                    odtcurrent_v72IR.warmt = tvIR.warmt;
214                    odtcurrent_v72IR.warmlatitude = tvIR.warmlatitude;
215                    odtcurrent_v72IR.warmlongitude = tvIR.warmlongitude;
216                    odtcurrent_v72IR.eyet = tvIR.eyet;
217                    odtcurrent_v72IR.cwcloudt = tvIR.cwcloudt;
218                    odtcurrent_v72IR.cwring = tvIR.cwring;
219    
220                    /*
221                     * Determine scene type Set scene type
222                     */
223    
224                    float[] oscen = StormAODTSceneType.aodtv72_calcscene(odtcurrent_v72IR,
225                                    areadata_v72);
226    
227                    odtcurrent_v72IR.cloudt = oscen[0];
228                    odtcurrent_v72IR.cloudt2 = oscen[1];
229                    odtcurrent_v72IR.eyestdv = oscen[2];
230                    odtcurrent_v72IR.cloudsymave = oscen[3];
231                    odtcurrent_v72IR.eyefft = (int) oscen[4];
232                    odtcurrent_v72IR.cloudfft = (int) oscen[5];
233                    // { alst, Aaveext, Estdveye, Aavesym, eyecnt, rngcnt};
234                    float[] oscen1 = StormAODTSceneType.aodtv72_classify(odtcurrent_v72IR,
235                                    rmwsizeman_v72, areadata_v72, osstr_v72, osearch_v72);
236    
237                    odtcurrent_v72IR.eyescene = (int) oscen1[0];
238                    odtcurrent_v72IR.cloudscene = (int) oscen1[1];
239                    odtcurrent_v72IR.eyesceneold = -1;
240                    odtcurrent_v72IR.cloudsceneold = -1;
241                    odtcurrent_v72IR.eyecdosize = oscen1[2];
242                    odtcurrent_v72IR.ringcb = (int) oscen1[3];
243                    odtcurrent_v72IR.ringcbval = (int) oscen1[4];
244                    odtcurrent_v72IR.ringcbvalmax = (int) oscen1[5];
245                    odtcurrent_v72IR.ringcblatmax = oscen1[6];
246                    odtcurrent_v72IR.ringcblonmax = oscen1[7];
247                    odtcurrent_v72IR.rmw = oscen1[8];
248    
249                    /*
250                     * Determine intensity
251                     */
252    
253                    iaodt = aodtv72_calcintensity(idomain);
254                    if (iaodt == 71) {
255                            throw new IllegalStateException("center location is over land");
256                    }
257    
258                    /*
259                     * Print out all diagnostic messages to screen
260                     */
261                    return odtcurrent_v72IR;
262            }
263    
264            public static class AODTResult {
265            }
266    
267            /**
268             * _more_
269             * 
270             * @param g2d
271             *            _more_
272             * @param cenlat
273             *            _more_
274             * @param cenlon
275             *            _more_
276             * @param numx
277             *            _more_
278             * @param numy
279             *            _more_
280             * 
281             * @return _more_
282             */
283    
284            GridUtil.Grid2D spatialSubset(GridUtil.Grid2D g2d, float cenlat,
285                            float cenlon, int numx, int numy) {
286                    float[][] lats = g2d.getlats();
287                    float[][] lons = g2d.getlons();
288                    float[][][] values = g2d.getvalues();
289                    float[][] slats = new float[numx][numy];
290                    float[][] slons = new float[numx][numy];
291                    float[][][] svalues = new float[1][numx][numy];
292    
293                    int ly = lats[0].length;
294                    int ly0 = ly / 2;
295                    int lx = lats.length;
296                    int lx0 = lx / 2;
297                    int ii = numx / 2, jj = numy / 2;
298    
299                    for (int j = 0; j < ly - 1; j++) {
300                            if (Float.isNaN(lats[lx0][j]))
301                                    continue;
302                            if ((lats[lx0][j] > cenlat) && (lats[lx0][j + 1] < cenlat)) {
303                                    jj = j;
304                            }
305                    }
306                    for (int i = 0; i < lx - 1; i++) {
307                            if (Float.isNaN(lons[i][ly0]))
308                                    continue;
309                            if ((lons[i][ly0] < cenlon) && (lons[i + 1][ly0] > cenlon)) {
310                                    ii = i;
311                            }
312                    }
313                    int startx = ii - (numx / 2 - 1);
314                    int starty = jj - (numy / 2 - 1);
315                    if (startx < 0)
316                            startx = 0;
317                    if (starty < 0)
318                            starty = 0;
319    
320                    for (int i = 0; i < numx; i++) {
321                            for (int j = 0; j < numy; j++) {
322                                    slats[i][j] = lats[i + startx][j + starty];
323                                    slons[i][j] = lons[i + startx][j + starty];
324                                    svalues[0][i][j] = values[0][i + startx][j + starty];
325                            }
326                    }
327    
328                    return new GridUtil.Grid2D(slats, slons, svalues);
329            }
330    
331            /**
332             * _more_
333             * 
334             * @param nx
335             *            _more_
336             * @param ny
337             *            _more_
338             * @param gv
339             *            _more_
340             * @param imsorc
341             *            _more_
342             * @param imtype
343             *            _more_
344             * 
345             * @return _more_
346             */
347            float[][] im_gvtota(int nx, int ny, float[][] gv, int imsorc, int imtype)
348    
349            /**
350             * im_gvtota
351             * 
352             * This subroutine converts GVAR counts to actual temperatures based on the
353             * current image set in IM_SIMG.
354             * 
355             * im_gvtota ( int *nvals, unsigned int *gv, float *ta, int *iret )
356             * 
357             * Input parameters: *nvals int Number of values to convert *gv int Array of
358             * GVAR count values
359             * 
360             * Output parameters: *ta float Array of actual temperatures *iret int
361             * Return value = -1 - could not open table = -2 - could not find match
362             * 
363             * 
364             * Log: D.W.Plummer/NCEP 02/03 D.W.Plummer/NCEP 06/03 Add coeff G for 2nd
365             * order poly conv T. Piper/SAIC 07/06 Added tmpdbl to eliminate warning
366             */
367            {
368                    int ii, ip, chan, found, ier;
369                    double Rad, Teff, tmpdbl;
370                    float[][] ta = new float[nx][ny];
371                    int iret;
372                    String fp = "/ucar/unidata/data/storm/ImgCoeffs.tbl";
373    
374                    iret = 0;
375    
376                    for (ii = 0; ii < nx; ii++) {
377                            for (int jj = 0; jj < ny; jj++) {
378                                    ta[ii][jj] = Float.NaN;
379                            }
380                    }
381    
382                    /*
383                     * Read in coefficient table if necessary.
384                     */
385                    String s = null;
386                    try {
387                            s = IOUtil.readContents(fp);
388                    } catch (Exception re) {
389                    }
390    
391                    int i = 0;
392                    StormAODTInfo.ImgCoeffs[] ImageConvInfo = new StormAODTInfo.ImgCoeffs[50];
393                    for (String line : StringUtil.split(s, "\n", true, true)) {
394                            if (line.startsWith("!")) {
395                                    continue;
396                            }
397                            List<String> stoks = StringUtil.split(line, " ", true, true);
398    
399                            ImageConvInfo[i] = new StormAODTInfo.ImgCoeffs(stoks);
400                            ;
401                            i++;
402                    }
403                    int nImgRecs = i;
404                    found = 0;
405                    ii = 0;
406                    while ((ii < nImgRecs) && (found == 0)) {
407    
408                            tmpdbl = (double) (ImageConvInfo[ii].chan - 1)
409                                            * (ImageConvInfo[ii].chan - 1);
410                            chan = G_NINT(tmpdbl);
411    
412                            if ((imsorc == ImageConvInfo[ii].sat_num) && (imtype == chan)) {
413                                    found = 1;
414                            } else {
415                                    ii++;
416                            }
417    
418                    }
419    
420                    if (found == 0) {
421                            iret = -2;
422                            return null;
423                    } else {
424    
425                            ip = ii;
426                            for (ii = 0; ii < nx; ii++) {
427                                    for (int jj = 0; jj < ny; jj++) {
428    
429                                            /*
430                                             * Convert GVAR count (gv) to Scene Radiance
431                                             */
432                                            Rad = ((double) gv[ii][jj] - ImageConvInfo[ip].scal_b) /
433                                            /* ------------------------------------- */
434                                            ImageConvInfo[ip].scal_m;
435    
436                                            Rad = Math.max(Rad, 0.0);
437    
438                                            /*
439                                             * Convert Scene Radiance to Effective Temperature
440                                             */
441                                            Teff = (c2 * ImageConvInfo[ip].conv_n)
442                                                            /
443                                                            /*
444                                                             * --------------------------------------------------
445                                                             * -----
446                                                             */
447                                                            (Math.log(1.0
448                                                                            + (c1 * Math.pow(ImageConvInfo[ip].conv_n,
449                                                                                            3.0)) / Rad));
450    
451                                            /*
452                                             * Convert Effective Temperature to Temperature
453                                             */
454                                            ta[ii][jj] = (float) (ImageConvInfo[ip].conv_a
455                                                            + ImageConvInfo[ip].conv_b * Teff + ImageConvInfo[ip].conv_g
456                                                            * Teff * Teff);
457                                    }
458    
459                            }
460    
461                    }
462    
463                    return ta;
464    
465            }
466    
467            /**
468             * _more_
469             * 
470             * @param x
471             *            _more_
472             * 
473             * @return _more_
474             */
475            public int G_NINT(double x) {
476                    return (((x) < 0.0F) ? ((((x) - (float) ((int) (x))) <= -.5f) ? (int) ((x) - .5f)
477                                    : (int) (x))
478                                    : ((((x) - (float) ((int) (x))) >= .5f) ? (int) ((x) + .5f)
479                                                    : (int) (x)));
480            }
481    
482            /**
483             * _more_
484             * 
485             * @param keyerM_v72
486             *            _more_
487             * @param areadata
488             *            _more_
489             * 
490             * @return _more_
491             */
492            StormAODTInfo.IRData aodtv72_seteyecloudtemp(int keyerM_v72,
493                            StormAODTInfo.DataGrid areadata)
494            /*
495             * Routine to search for, identify, and set the eye and cloud temperature
496             * values for the AODT library. Temperatures are set within AODT library.
497             * Inputs : none Outputs: none Return : -51 : eye, CWcloud, or warmest
498             * temperature <-100C or >+40C 0 : o.k.
499             */
500            {
501                    StormAODTInfo.IRData ird = StormAODTSceneType.aodtv72_gettemps(
502                                    keyerM_v72, areadata);
503                    if (ird == null) {
504                            throw new IllegalStateException(
505                                            "eye, CWcloud, or warmest temperature <-100C or >+40C");
506                    }
507    
508                    return ird;
509    
510                    // return iok;
511            }
512    
513            /**
514             * _more_
515             * 
516             * @param temps
517             *            _more_
518             * @param lats
519             *            _more_
520             * @param lons
521             *            _more_
522             * @param numx
523             *            _more_
524             * @param numy
525             *            _more_
526             * 
527             * @return _more_
528             */
529            public int aodtv72_loadIRimage(float[][] temps, float[][] lats,
530                            float[][] lons, int numx, int numy)
531            /*
532             * Subroutine to load IR image data grid values (temperatures and positions)
533             * into data structure for AODT library Inputs : temperature, latitude, and
534             * longitude arrays centered on storm position location along with number of
535             * columns (x) and rows (y) in grid Outputs: none (areadata_v72 structure
536             * passed via global variable) Return : 0 : o.k.
537             */
538            {
539                    StormAODTInfo sinfo = new StormAODTInfo();
540                    /* allocate space for data */
541                    areadata_v72 = new StormAODTInfo.DataGrid(temps, lats, lons, numx, numy);
542                    return 0;
543            }
544    
545            /**
546             * _more_
547             * 
548             * @param indomain
549             *            _more_
550             * 
551             * @return _more_
552             */
553            int aodtv72_setdomain(int indomain)
554            /*
555             * set current ocean domain variable within AODT library memory Inputs :
556             * domain flag value from input Outputs: none Return : -81 : error
557             * determining storm basin
558             */
559            {
560                    int domain;
561                    float xlon;
562    
563                    /* obtain current storm center longitude */
564                    xlon = odtcurrent_v72IR.longitude;
565                    if ((xlon < -180.0) || (xlon > 180.0)) {
566                            return -81;
567                    }
568    
569                    ixdomain_v72 = indomain;
570                    /* determine oceanic domain */
571                    if (indomain == 0) {
572                            /* automatically determined storm basin */
573                            if (xlon >= 0.0) {
574                                    domain = 0; /* atlantic and east pacific to 180W/dateline */
575                            } else {
576                                    domain = 1; /* west pacific and other regions */
577                            }
578                    } else {
579                            /* manually determined storm basin */
580                            domain = indomain - 1;
581                    }
582    
583                    /* assign ocean domain flag value to AODT library variable */
584                    idomain_v72 = domain;
585    
586                    return 0;
587            }
588    
589            /**
590             * _more_
591             * 
592             * @param ilat
593             *            _more_
594             * @param ilon
595             *            _more_
596             * @param ipos
597             *            _more_
598             * 
599             * @return _more_
600             */
601            int aodtv72_setlocation(float ilat, float ilon, int ipos)
602            /*
603             * set current storm center location within from AODT library memory Inputs
604             * : AODT library current storm center latitude and longitude values and
605             * location positioning method : 1-forecast interpolation 2-laplacian
606             * technique 3-warm spot 4-extrapolation Outputs: none Return : -21 :
607             * invalid storm center position 21 : user selected storm center position 22
608             * : auto selected storm center position
609             */
610            {
611                    int iret;
612    
613                    /* assign current storm center latitude value to AODT library variable */
614                    odtcurrent_v72IR.latitude = ilat;
615                    /* assign current storm center longitude value to AODT library variable */
616                    odtcurrent_v72IR.longitude = ilon;
617                    /* assign current storm center positioning flag to AODT library variable */
618                    odtcurrent_v72IR.autopos = ipos;
619                    if ((odtcurrent_v72IR.longitude < -180.)
620                                    || (odtcurrent_v72IR.longitude > 180.)) {
621                            iret = -21;
622                    }
623                    if ((odtcurrent_v72IR.latitude < -90.)
624                                    || (odtcurrent_v72IR.latitude > 90.)) {
625                            iret = -21;
626                    }
627    
628                    iret = 21; /* user selected image location */
629                    if (ipos >= 1) {
630                            iret = 22;
631                    }
632    
633                    return iret;
634            }
635    
636            /**
637             * _more_
638             * 
639             * @param datetime
640             *            _more_
641             * @param sat
642             *            _more_
643             * 
644             * @return _more_
645             */
646            int aodtv72_setIRimageinfo(double datetime, int sat)
647            /*
648             * set IR image date/time within AODT library memory Inputs : AODT library
649             * IR image date/time/satellite information Outputs: none Return : 0 : o.k.
650             */
651            {
652                    /* assign IR image date to AODT library variable */
653                    odtcurrent_v72IR.date = datetime;
654    
655                    /* assign IR image satellite type to AODT library variable */
656                    odtcurrent_v72IR.sattype = sat;
657    
658                    return 0;
659            }
660    
661            /**
662             * _more_
663             * 
664             * @param idomain
665             *            _more_
666             * 
667             * @return _more_
668             */
669            public int aodtv72_calcintensity(int idomain)
670            /*
671             * Compute intensity values CI, Final T#, and Raw T#. Inputs : global
672             * structure odtcurrent_v72 containing current analysis Outputs : none
673             * Return : 71 : storm is over land 0 : o.k.
674             */
675            {
676                    int iok;
677                    int iret;
678                    int strength;
679    
680                    if ((odtcurrent_v72IR.land == 1)) {
681                            iok = aodtv72_initcurrent(true, odtcurrent_v72IR);
682                            iret = 71;
683                    } else {
684                            /* calculate current Raw T# value */
685                            odtcurrent_v72IR.Traw = aodtv72_Tnoraw(odtcurrent_v72IR, idomain);
686                            odtcurrent_v72IR.TrawO = odtcurrent_v72IR.Traw;
687                            /* check for spot analysis or full analysis using history file */
688                            /* if(hfile_v72==(char *)NULL) { */
689                            if (true) {
690                                    /* perform spot analysis (only Traw) */
691                                    odtcurrent_v72IR.Tfinal = odtcurrent_v72IR.Traw;
692                                    odtcurrent_v72IR.Tfinal3 = odtcurrent_v72IR.Traw;
693                                    odtcurrent_v72IR.CI = odtcurrent_v72IR.Traw;
694                                    odtcurrent_v72IR.CIadjp = aodtv72_latbias(odtcurrent_v72IR.CI,
695                                                    odtcurrent_v72IR.latitude, odtcurrent_v72IR.longitude,
696                                                    odtcurrent_v72IR);
697                                    /*
698                                     * printf("%f %f %f   %f\n",odtcurrent_v72IR.CI,odtcurrent_v72IR.
699                                     * latitude
700                                     * ,odtcurrent_v72->IR.longitude,odtcurrent_v72->IR.CIadjp);
701                                     */
702                                    odtcurrent_v72IR.rule9 = 0;
703                                    /* odtcurrent_v72->IR.TIEraw=aodtv72_TIEmodel(); */
704                                    /* odtcurrent_v72->IR.TIEavg=odtcurrent_v72->IR.TIEraw; */
705                                    /* odtcurrent_v72->IR.TIEflag=aodtv72_tieflag(); */
706                            } else {
707                            }
708    
709                            iret = 0;
710                    }
711    
712                    return iret;
713            }
714    
715            /**
716             * _more_
717             * 
718             * @param initval
719             *            _more_
720             * @param latitude
721             *            _more_
722             * @param longitude
723             *            _more_
724             * @param odtcurrent_v72IR
725             *            _more_
726             * 
727             * @return _more_
728             */
729            float aodtv72_latbias(float initval, float latitude, float longitude,
730                            StormAODTInfo.IRData odtcurrent_v72IR)
731            /*
732             * Apply Latitude Bias Adjustment to CI value Inputs : initval - initial CI
733             * value latitude - current latitude of storm Outputs : adjusted MSLP value
734             * as return value
735             */
736            {
737                    float initvalp;
738    
739                    float initvalpx;
740                    float value; /* lat bias adjustement amount (0.00-1.00) */
741                    int sceneflag; /*
742                                                     * contains lat bias adjustment flag 0=no adjustment
743                                                     * 1=intermediate adjustment (6 hours) 2=full adjustment
744                                                     */
745    
746                    sceneflag = aodtv72_scenesearch(0); /*
747                                                                                             * 0 means search for EIR based
748                                                                                             * parameters... cdo, etc
749                                                                                             */
750                    value = 1.0f; /* this value should be return from scenesearch() */
751                    /* printf("sceneflag=%d  value=%f\n",sceneflag,value); */
752                    odtcurrent_v72IR.LBflag = sceneflag;
753                    /* initvalp=aodtv72_getpwval(0,initval); TLO */
754                    initvalp = 0.0f;
755                    if (sceneflag >= 2) {
756                            /* EIR scene */
757                            if ((latitude >= 0.0)
758                                            && ((longitude >= -100.0) && (longitude <= -40.0))) {
759                                    /* do not make adjustment in N Indian Ocean */
760                                    return initvalp;
761                            }
762                            /* apply bias adjustment to pressure */
763                            /* initvalp=-1.0*value*(-20.60822+(0.88463*A_ABS(latitude))); */
764                            initvalp = value * (7.325f - (0.302f * Math.abs(latitude)));
765                    }
766    
767                    return initvalp;
768            }
769    
770            /**
771             * _more_
772             * 
773             * @param type
774             *            _more_
775             * 
776             * @return _more_
777             */
778            int aodtv72_scenesearch(int type) {
779                    int curflag = 1, flag, eirflag;
780                    float eirvalue, civalue, ciadjvalue, latitude;
781                    double curtime, xtime, curtimem6, mergetimefirst, mergetimelast, firsttime = -9999.0;
782    
783                    float pvalue;
784    
785                    /*
786                     * if(((odthistoryfirst_v72==0)&&(ostartstr_v72==TRUE))&&(hfile_v72!=(char
787                     * *)NULL)) {
788                     */
789    
790                    if (true) {
791                            flag = 2;
792                            pvalue = 1.0f;
793                            return flag;
794                    }
795    
796                    return flag;
797            }
798    
799            /**
800             * _more_
801             * 
802             * @param redo
803             *            _more_
804             * @param odtcurrent_v72IR
805             *            _more_
806             * 
807             * @return _more_
808             */
809            int aodtv72_initcurrent(boolean redo, StormAODTInfo.IRData odtcurrent_v72IR)
810            /*
811             * initialize odtcurrent_v72 array or reset values for land interaction
812             * situations
813             */
814            {
815                    int b, bb;
816                    // char comm[50]="\0";
817    
818                    if (!redo) {
819                            // odtcurrent_v72=(struct odtdata *)malloc(sizeof(struct odtdata));
820                            odtcurrent_v72IR.latitude = 999.99f;
821                            odtcurrent_v72IR.longitude = 999.99f;
822                            odtcurrent_v72IR.land = 0;
823                            odtcurrent_v72IR.autopos = 0;
824                            // strcpy(odtcurrent_v72IR.comment,comm);
825                            // diagnostics_v72=(char *)calloc((size_t)50000,sizeof(char));
826                            // hfile_v72=(char *)calloc((size_t)200,sizeof(char));
827                            // fixfile_v72=(char *)calloc((size_t)200,sizeof(char));
828    
829                            // b=sizeof(float);
830                            // bb=sizeof(double);
831                            // fcstlat_v72=(float *)calloc((size_t)5,b);
832                            // fcstlon_v72=(float *)calloc((size_t)5,b);
833                            // fcsttime_v72=(double *)calloc((size_t)5,bb);
834                    }
835    
836                    odtcurrent_v72IR.Traw = 0.0f;
837                    odtcurrent_v72IR.TrawO = 0.0f;
838                    odtcurrent_v72IR.Tfinal = 0.0f;
839                    odtcurrent_v72IR.Tfinal3 = 0.0f;
840                    odtcurrent_v72IR.CI = 0.0f;
841                    odtcurrent_v72IR.eyet = 99.99f;
842                    odtcurrent_v72IR.warmt = 99.99f;
843                    odtcurrent_v72IR.cloudt = 99.99f;
844                    odtcurrent_v72IR.cloudt2 = 99.99f;
845                    odtcurrent_v72IR.cwcloudt = 99.99f;
846                    odtcurrent_v72IR.warmlatitude = 999.99f;
847                    odtcurrent_v72IR.warmlongitude = 999.99f;
848                    odtcurrent_v72IR.eyecdosize = 0.0f;
849                    odtcurrent_v72IR.eyestdv = 0.0f;
850                    odtcurrent_v72IR.cloudsymave = 0.0f;
851                    odtcurrent_v72IR.eyescene = 0;
852                    odtcurrent_v72IR.cloudscene = 0;
853                    odtcurrent_v72IR.eyesceneold = -1;
854                    odtcurrent_v72IR.cloudsceneold = -1;
855                    odtcurrent_v72IR.rule9 = 0;
856                    odtcurrent_v72IR.rule8 = 0;
857                    odtcurrent_v72IR.LBflag = 0;
858                    odtcurrent_v72IR.rapiddiss = 0;
859                    odtcurrent_v72IR.eyefft = 0;
860                    odtcurrent_v72IR.cloudfft = 0;
861                    odtcurrent_v72IR.cwring = 0;
862                    odtcurrent_v72IR.ringcb = 0;
863                    odtcurrent_v72IR.ringcbval = 0;
864                    odtcurrent_v72IR.ringcbvalmax = 0;
865                    odtcurrent_v72IR.CIadjp = 0.0f;
866                    odtcurrent_v72IR.rmw = -99.9f;
867                    /* odtcurrent_v72->IR.TIEflag=0; */
868                    /* odtcurrent_v72->IR.TIEraw=0.0; */
869                    /* odtcurrent_v72->IR.TIEavg=0.0; */
870                    /* odtcurrent_v72->IR.sst=-99.9; */
871                    // if(!redo) odtcurrent_v72->nextrec=NULL; /* added by CDB */
872    
873                    return 0;
874            }
875    
876            /**
877             * Compute initial Raw T-Number value using original Dvorak rules
878             * 
879             * @param odtcurrent
880             * @param idomain_v72
881             * @return return value is Raw T#
882             */
883    
884            float aodtv72_Tnoraw(StormAODTInfo.IRData odtcurrent, int idomain_v72)
885            /*
886             * Compute initial Raw T-Number value using original Dvorak rules Inputs :
887             * global structure odtcurrent_v72 containing current analysis Outputs :
888             * return value is Raw T#
889             * 
890             * ODT SCENE/TEMPERATURE TABLE BD | WMG OW DG MG LG B W CMG CDG | TEMP |30.0
891             * 0.0 -30.0 -42.0 -54.0 -64.0 -70.0 -76.0 -80.0+|
892             * ---------------------------------------------------------------| Atl EYE
893             * | 3.5 4.0 4.5 4.5 5.0 5.5 6.0 6.5 7.0 | EMBC | 3.5 3.5 4.0 4.0 4.5 4.5
894             * 5.0 5.0 5.0 | CDO | 3.0 3.0 3.5 4.0 4.5 4.5 4.5 5.0 5.0 |
895             * ---------------------------------------------------------------| Pac EYE
896             * | 4.0 4.0 4.0 4.5 4.5 5.0 5.5 6.0 6.5 | EMBC | 3.5 3.5 4.0 4.0 4.5 4.5
897             * 5.0 5.0 5.0 | CDO | 3.0 3.5 3.5 4.0 4.5 4.5 4.5 4.5 5.0 |
898             * ---------------------------------------------------------------| Cat diff
899             * | 0 1 2 3 4 5 6 7 8 | add | 0.0 0.0 0.0 0.0 0.0-->0.5 0.5-->1.0 1.5 |
900             * (old) add |-0.5 -0.5 0.0 0.0-->0.5 0.5 0.5-->1.0 1.0 | (new)
901             * ---------------------------------------------------------------|
902             */
903            {
904    
905                    double[][] eno = {
906                                    { 1.00, 2.00, 3.25, 4.00, 4.75, 5.50, 5.90, 6.50, 7.00, 7.50,
907                                                    8.00 }, /* original plus adjusted > CDG+ */
908                                    { 1.50, 2.25, 3.30, 3.85, 4.50, 5.00, 5.40, 5.75, 6.25, 6.50,
909                                                    7.00 } }; /* adjusted based */
910                    double[][] cdo = {
911                                    { 2.00, 2.40, 3.25, 3.50, 3.75, 4.00, 4.10, 4.20, 4.30, 4.40,
912                                                    4.70 },
913                                    { 2.05, 2.40, 3.00, 3.20, 3.40, 3.55, 3.65, 3.75, 3.80, 3.90,
914                                                    4.10 } };
915                    double[] curbnd = { 1.0, 1.5, 2.5, 3.0, 3.5, 4.0, 4.5 };
916                    double[] shrdst = { 0.0, 35.0, 50.0, 80.0, 110.0, 140.0 };
917                    double[] shrcat = { 3.5, 3.0, 2.5, 2.25, 2.0, 1.5 };
918    
919                    double[][] diffchk = {
920                                    { 0.0, 0.5, 1.2, 1.7, 2.2, 2.7, 0.0, 0.0, 0.1, 0.5 }, /*
921                                                                                                                                             * shear
922                                                                                                                                             * scene
923                                                                                                                                             * types...
924                                                                                                                                             * original
925                                                                                                                                             * Rule 8
926                                                                                                                                             * rules
927                                                                                                                                             */
928                                    { 0.0, 0.5, 1.7, 2.2, 2.7, 3.2, 0.0, 0.0, 0.1, 0.5 }, /*
929                                                                                                                                             * eye scene
930                                                                                                                                             * types...
931                                                                                                                                             * add 0.5
932                                                                                                                                             * to Rule 8
933                                                                                                                                             * rules
934                                                                                                                                             */
935                                    { 0.0, 0.5, 0.7, 1.2, 1.7, 2.2, 0.0, 0.0, 0.1, 0.5 } }; /*
936                                                                                                                                                     * other
937                                                                                                                                                     * scene
938                                                                                                                                                     * types
939                                                                                                                                                     * ...
940                                                                                                                                                     * subtract
941                                                                                                                                                     * 0.5
942                                                                                                                                                     * from
943                                                                                                                                                     * Rule
944                                                                                                                                                     * 8
945                                                                                                                                                     * rules
946                                                                                                                                                     */
947                    double eyeadjfacEYE[] = { 0.011, 0.015 }; /*
948                                                                                                     * modified wpac value to be
949                                                                                                     * closer to atlantic
950                                                                                                     */
951                    double symadjfacEYE[] = { -0.015, -0.015 };
952                    double dgraysizefacCLD[] = { 0.002, 0.001 };
953                    double symadjfacCLD[] = { -0.030, -0.015 };
954    
955                    int diffchkcat;
956                    int ixx, cloudcat, eyecat, diffcat, rp, xrp, rb;
957                    float incval, lastci, lasttno, lastr9, lastraw;
958                    float xpart, xparteye, xaddtno, eyeadj, spart, ddvor, dvorchart, ciadj;
959                    float sdist, cloudtemp, eyetemp, fftcloud;
960                    float t1val, t6val, t12val, t18val, t24val, delt1, delt6, delt12, delt18, delt24;
961                    float t1valraw, t1valrawx, txvalmin, txvalmax;
962                    double curtime, xtime, firsttime, firstlandtime;
963                    double ttime1, ttime6, ttime12, ttime18, ttime24, t1valrawxtime;
964                    StormAODTInfo.IRData odthistory, prevrec;
965                    boolean oceancheck, adjustshear, firstland;
966                    boolean t1found = false, t6found = false, t12found = false, t18found = false, t24found = false;
967                    boolean first6hrs = false;
968                    float symadj, dgraysizeadj, deltaT;
969    
970                    cloudtemp = odtcurrent.cloudt;
971                    eyetemp = odtcurrent.eyet;
972                    cloudcat = 0;
973                    eyecat = 0;
974                    lastci = 4.0f;
975                    xpart = 0.0f;
976    
977                    for (ixx = 0; ixx < 10; ixx++) {
978                            /* compute cloud category */
979                            if ((cloudtemp <= StormAODTInfo.ebd_v72[ixx])
980                                            && (cloudtemp > StormAODTInfo.ebd_v72[ixx + 1])) {
981                                    cloudcat = ixx;
982                                    xpart = (float) (cloudtemp - StormAODTInfo.ebd_v72[cloudcat])
983                                                    / (float) (StormAODTInfo.ebd_v72[cloudcat + 1] - StormAODTInfo.ebd_v72[cloudcat]);
984                            }
985                            /* compute eye category for eye adjustment */
986                            if ((eyetemp <= StormAODTInfo.ebd_v72[ixx])
987                                            && (eyetemp > StormAODTInfo.ebd_v72[ixx + 1])) {
988                                    eyecat = ixx;
989                            }
990                            /* eyetemp=Math.min(0.0,eyetemp); */
991                    }
992                    if (odtcurrent.eyescene == 1) {
993                            /* for pinhole eye, determine what storm should be seeing */
994                            /*
995                             * eyetemp=pinhole(odtcurrent_v72->IR.latitude,odtcurrent_v72->IR.longitude
996                             * ,eyetemp);
997                             */
998                            /*
999                             * eyetemp=(9.0-eyetemp)/2.0; / this matches DT used at NHC (jack
1000                             * beven)
1001                             */
1002                            eyetemp = (float) (eyetemp - 9.0) / 2.0f; /*
1003                                                                                                             * between +9C (beven) and
1004                                                                                                             * measured eye temp (turk)
1005                                                                                                             */
1006                            odtcurrent.eyet = eyetemp;
1007                    }
1008    
1009                    /* category difference between eye and cloud region */
1010                    diffcat = Math.max(0, cloudcat - eyecat);
1011    
1012                    /* if scenetype is EYE */
1013                    rp = odtcurrent.ringcbval;
1014                    rb = odtcurrent.ringcb;
1015                    fftcloud = odtcurrent.cloudfft;
1016    
1017                    if (odtcurrent.cloudscene == 3) {
1018                            /* CURVED BAND */
1019                            rp = Math.min(30, rp + 1); /* added 1 for testing */
1020                            xrp = rp / 5;
1021                            incval = 0.1f;
1022                            if (xrp == 1) {
1023                                    incval = 0.2f;
1024                            }
1025                            ddvor = (float) curbnd[xrp];
1026                            xaddtno = incval * (float) (rp - (xrp * 5));
1027                            /*
1028                             * printf("rp=%d  xrp=%d  rb=%d  ddvor=%f  xaddtno=%f\n",rp,xrp,rb,ddvor
1029                             * ,xaddtno);
1030                             */
1031                            ddvor = ddvor + xaddtno;
1032                            if (rb == 5) {
1033                                    ddvor = Math.min(4.0f, ddvor + 0.5f);
1034                            }
1035                            if (rb == 6) {
1036                                    ddvor = Math.min(4.5f, ddvor + 1.0f);
1037                            }
1038                            diffchkcat = 2; /* added for test - non-eye/shear cases */
1039                    } else if (odtcurrent.cloudscene == 4) {
1040                            /* POSSIBLE SHEAR -- new definition from NHC */
1041                            ixx = 0;
1042                            ddvor = 1.0f;
1043                            sdist = odtcurrent.eyecdosize; /* shear distance */
1044                            while (ixx < 5) {
1045                                    if ((sdist >= shrdst[ixx]) && (sdist < shrdst[ixx + 1])) {
1046                                            spart = (float) ((sdist - shrdst[ixx]) / (shrdst[ixx + 1] - shrdst[ixx]));
1047                                            xaddtno = (float) ((spart * (shrcat[ixx + 1] - shrcat[ixx])));
1048                                            ddvor = (float) (shrcat[ixx] + xaddtno);
1049                                            ixx = 5;
1050                                    } else {
1051                                            ixx++;
1052                                    }
1053                            }
1054                            diffchkcat = 0; /* added for test - shear cases */
1055                    } else {
1056                            /* EYE or NO EYE */
1057                            if (odtcurrent.eyescene <= 2) {
1058                                    /* EYE */
1059                                    xaddtno = (float) (xpart * (eno[idomain_v72][cloudcat + 1] - eno[idomain_v72][cloudcat]));
1060                                    /*
1061                                     * cloud category must be white (-70C) or below for full
1062                                     * adjustment; value will be merged in starting at black (-64C)
1063                                     * / if(cloudcat<5) { / gray shades / xparteye=0.00; } else
1064                                     * if(cloudcat==5) { / black / xparteye=xpart; } else { / white
1065                                     * and colder / xparteye=1.00; }
1066                                     */
1067                                    eyeadj = (float) eyeadjfacEYE[idomain_v72]
1068                                                    * (eyetemp - cloudtemp);
1069                                    /* symadj=-0.02*(odtcurrent_v72->IR.cloudsymave); */
1070                                    symadj = (float) symadjfacEYE[idomain_v72]
1071                                                    * (odtcurrent.cloudsymave);
1072                                    /*
1073                                     * printf("EYE : cloudsymave=%f  symadj=%f\n",odtcurrent_v72->IR.
1074                                     * cloudsymave,symadj);
1075                                     */
1076                                    ddvor = (float) eno[idomain_v72][cloudcat] + xaddtno + eyeadj
1077                                                    + symadj;
1078                                    /*
1079                                     * printf("EYE : xaddtno=%f  eyeadj=%f  symadj=%f   ddvor=%f\n",xaddtno
1080                                     * ,eyeadj,symadj,ddvor);
1081                                     */
1082                                    ddvor = Math.min(ddvor, 9.0f);
1083                                    /* printf("ddvor=%f\n",ddvor); */
1084                                    if (odtcurrent.eyescene == 2) {
1085                                            ddvor = Math.min(ddvor - 0.5f, 6.5f); /* LARGE EYE adjustment */
1086                                    }
1087                                    /*
1088                                     * if(odtcurrent_v72->IR.eyescene==3)
1089                                     * ddvor=Math.min(ddvor-0.5,6.0); / LARGE RAGGED EYE adjustment
1090                                     */
1091                                    diffchkcat = 1; /* added for test - eye cases */
1092                                    /* printf("ddvor=%f\n",ddvor); */
1093                            } else {
1094                                    /* NO EYE */
1095                                    /* CDO */
1096                                    xaddtno = (float) (xpart * (cdo[idomain_v72][cloudcat + 1] - cdo[idomain_v72][cloudcat]));
1097                                    /* dgraysizeadj=0.002*odtcurrent_v72->IR.eyecdosize; */
1098                                    dgraysizeadj = (float) dgraysizefacCLD[idomain_v72]
1099                                                    * odtcurrent.eyecdosize;
1100                                    /*
1101                                     * printf("CDO : dgraysize=%f  symadj=%f\n",odtcurrent_v72->IR.eyecdosize
1102                                     * ,dgraysizeadj);
1103                                     */
1104                                    /* symadj=-0.03*(odtcurrent_v72->IR.cloudsymave); */
1105                                    symadj = (float) symadjfacCLD[idomain_v72]
1106                                                    * (odtcurrent.cloudsymave);
1107                                    /*
1108                                     * printf("CDO : cloudsymave=%f  symadj=%f\n",odtcurrent_v72->IR.
1109                                     * cloudsymave,symadj);
1110                                     */
1111                                    ddvor = (float) cdo[idomain_v72][cloudcat] + xaddtno
1112                                                    + dgraysizeadj + symadj;
1113                                    ddvor = ddvor - 0.1f; /* bias adjustment */
1114                                    /*
1115                                     * printf("CDO : xaddtno=%f dgraysizeadj=%f  symadj=%f   ddvor=%f\n"
1116                                     * ,xaddtno,dgraysizeadj,symadj,ddvor);
1117                                     */
1118                                    ciadj = 0.0f;
1119                                    if (odtcurrent.cloudscene == 0) { /* CDO */
1120                                            if (lastci >= 4.5) {
1121                                                    ciadj = Math.max(0.0f, Math.min(1.0f, lastci - 4.5f));
1122                                            }
1123                                            if (lastci <= 3.0) {
1124                                                    ciadj = Math.min(0.0f, Math.max(-1.0f, lastci - 3.0f));
1125                                            }
1126                                            /* printf("CDO : lastci=%f   xaddtno=%f\n",lastci,ciadj); */
1127                                            ddvor = ddvor + ciadj;
1128                                    }
1129                                    if (odtcurrent.cloudscene == 1) { /* EMBEDDED CENTER */
1130                                            ciadj = Math.max(0.0f, Math.min(1.5f, lastci - 4.0f));
1131                                            /* printf("EMBC : lastci=%f   xaddtno=%f\n",lastci,ciadj); */
1132                                            ddvor = ddvor + ciadj; /* changed from 0.5 */
1133                                    }
1134                                    if (odtcurrent.cloudscene == 2) { /* IRREGULAR CDO (PT=3.5) */
1135                                            ddvor = ddvor + 0.3f; /* additional IrrCDO bias adjustment */
1136                                            ddvor = Math.min(3.5f, Math.max(2.5f, ddvor));
1137                                    }
1138                                    diffchkcat = 2; /* added for test - non-eye/shear cases */
1139                            }
1140                    }
1141    
1142                    dvorchart = ((float) (int) (ddvor * 10.0f)) / 10.0f;
1143                    // odtcurrent_v72IR.TrawO=dvorchart;
1144    
1145                    return dvorchart;
1146    
1147            }
1148    
1149    }