001/*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2017
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
029package edu.wisc.ssec.mcidasv.adt;
030
031import static java.lang.Math.PI;
032import static java.lang.Math.abs;
033import static java.lang.Math.cos;
034import static java.lang.Math.log;
035import static java.lang.Math.max;
036import static java.lang.Math.sqrt;
037import static java.lang.Math.pow;
038
039import java.io.IOException;
040
041import javax.swing.JOptionPane;
042
043class Remap {
044   /** Block size (bytes) input file */
045   int in_bfw;
046   
047   /** Block size (bytes) output file */
048   int out_bfw;
049   
050   /** Number of splines/line */
051   int nspl;
052   
053   /** Number of splines/elem */
054   int nspe;
055   
056   /** Source blocksize */
057   int slb;
058   
059   /** Dest blocksize */
060   int dlb;
061   
062   /** Number of corners in line */
063   int ncl;
064   
065   /** Number of corners in elem */
066   int nce;
067}
068
069/** TIFF header */
070class TiffHeader {
071   /** Byte order */
072   int order;
073   
074   /** Version */
075   int version;
076   
077   /** Pointer */
078   int point;
079}
080
081/**
082 * TIFF directory record.
083 */
084class TiffRecord {
085   /** TIFF tag */
086   int tag;
087   
088   /** Data type */
089   int type;
090   
091   /** Length */
092   int length;
093   
094   /** Pointer or value */
095   int voff;
096}
097
098class DisVars {
099   /** Output number of lines */
100   double xrectl;
101   
102   /** Output number of elems */
103   double xrecte;
104}
105
106class TiffVars {
107        int nbits;
108        int photo;
109        int unit;
110        int in_lines;
111        int in_elems;
112        int out_lines;
113        int out_elems;
114}
115
116public class Auto {
117
118   private static int[][] MoatMaskFlagField = new int[200][200];
119   
120   private static int[][] BlackWhiteFieldArray = new int[200][200];
121   
122   private static double[][] IRData_Remap_Latitude = new double[200][200];
123   private static double[][] IRData_Remap_Longitude = new double[200][200];
124   
125   private static double[][] IRData_Remap_Temperature = new double[200][200];
126   
127   private static int IRData_Remap_NumberRows;
128   
129   private static int IRData_Remap_NumberColumns;
130   
131   private static double[][] NSTempGradientArray = new double[200][200];
132   
133   private static double[][] EWTempGradientArray = new double[200][200];
134   
135   private static double[][] SpiralCenterAnalysisField = new double[50000][3];
136   
137   private static double[][] RingScoreAnalysisField = new double[750000][3];
138   
139   private static int[][] CircleFilterRowArray = new int[100][2000];
140   
141   private static int[][] CircleFilterColumnArray = new int[100][2000];
142   
143   /* I don't know this size for sure */
144   /** Array containing line coordinates */
145   private static double[] LineCoordinateArray = new double[5000];
146   
147   /* I don't know this size for sure */
148   /** Array containing element coordinates */
149   private static double[] ElementCoordinateArray = new double[5000];
150
151   private static Remap remap_vars = new Remap();
152   
153   private static TiffVars tiff_vars = new TiffVars();
154   
155   private static double RING_WIDTH = 4.0;
156   
157   /** Minimum block size (bytes) for domap */
158   private static int MINBFW = 500000;
159   
160   /** Minimum block size (lines) */
161   private static int MINBLKSIZ = 10;
162
163   public Auto() {
164     IRData_Remap_NumberRows = 0;
165     IRData_Remap_NumberColumns = 0;
166   }
167   
168   /**
169    * Determine storm position at time CurrentTime using NHC/JTWC
170    * forecast discussion products.
171    *
172    * Time and location information
173    * from these products are then interpolated to time in question
174    * to derive a estimated storm position.  If position estimation
175    * cannot be calculated, a lat/lon position of -99.5/-999.5 will
176    * be returned.
177    * Inputs  : None
178    * Outputs : Latitude_Return            - estimated latitude position
179    *           Longitude_Return           - estimated longitude position
180    *           PositioningMethodID_Return - method used to derive storm location
181    * Return  : -43 : Error w/ forecast file open and BAD extrapolation
182    *           -44 : Invalid forecast file and BAD extrapolation
183    *           -45 : Error w/ forecast file read and BAD extrapolation
184    *           -46 : Error w/ forecast interpolation and BAD extrapolation
185    *            42 : GOOD INTERPOLATION
186    *            43 : Error w/ forecast file open but GOOD EXTRAPOLATION
187    *            44 : Invalid forecast file but GOOD extrapolation
188    *            45 : Error w/ forecast file read but GOOD extrapolation
189    *            46 : Error w/ forecast interpolation but GOOD EXTRAPOLATION
190    *             0 : Subroutine Error
191    */
192   public double[] AutoMode1(String ForecastFile, int ForecastFileType) throws IOException {
193
194
195      int RetID = 0;
196      int PositioningMethodID = 0;
197      int InterpolatedReturnFlag = 0;
198      double InterpolatedLatitude = -99.99;
199      double InterpolatedLongitude = -99.99;
200      double InterpolatedIntensity = -99.99;
201      boolean UseExtrapPositionTF = false;
202
203      /* Read Forecast File */
204      double ThresholdTime = 24.0;
205      double[] ForecastFileOutput = Forecasts.ReadForecasts(ForecastFile,ForecastFileType,ThresholdTime);
206      
207      if (ForecastFileOutput == null) {
208          JOptionPane.showMessageDialog(null, "Unable to process Forecast File, please ensure a valid file is present.");
209          return null;
210      }
211      
212      InterpolatedReturnFlag = (int)ForecastFileOutput[0];
213      InterpolatedLatitude = ForecastFileOutput[1];
214      InterpolatedLongitude = ForecastFileOutput[2];
215      InterpolatedIntensity = ForecastFileOutput[3];
216      /* History.IRCurrentRecord.latitude = ForecastLatitude; */
217      /* History.IRCurrentRecord.longitude = ForecastLongitude; */
218
219      /* System.out.printf("InterpolatedReturnFlag=%d\n",InterpolatedReturnFlag); */
220      if(InterpolatedReturnFlag!=0) {
221         PositioningMethodID = 6;
222         UseExtrapPositionTF = true;
223         /*
224          * -1 will give 45=invalid interpretation
225          * -2 will give 44=invalid file type
226          * -4 will give 42=error closing file
227          * -5 will give 41=error opening file
228          */
229         RetID = 46+InterpolatedReturnFlag;
230      } else {
231         /* Good interpolated values... use 'em */
232         PositioningMethodID = 1;
233         UseExtrapPositionTF = false;
234         RetID = 43;
235      }
236
237      /*
238       * try to extrapolate storm location from previous
239       *  storm locations in history file
240       */
241      if(UseExtrapPositionTF) {
242         /* call slopecal to get y-intercept values for lat/lon values */
243         InterpolatedLatitude = Functions.adt_slopecal(12.0,3);
244         InterpolatedLongitude = Functions.adt_slopecal(12.0,4);
245         if((abs(InterpolatedLatitude)>90.0)||
246            (abs(InterpolatedLongitude)>180.0)) {
247            /* invalid interp and extrap... negative error code returns */
248            PositioningMethodID = 0;
249            RetID = -1*RetID;
250         } else {
251            /* invalid interp but good extrap... positive error code returns */
252            PositioningMethodID = 6;
253            RetID = 46;
254         }
255      }
256
257      return new double[] {
258         (double)RetID,
259         InterpolatedLatitude,
260         InterpolatedLongitude,
261         InterpolatedIntensity,
262         (double)PositioningMethodID
263      };
264   }
265   
266   /**
267    * Additional automatic positioning of storm center location using
268    * official forecasts from NHC or JTWC as input.
269    *
270    * Storm location will be estimated using spiral fitting and ring fitting
271    * routines derived by Tony Wimmers in his MatLab routines (which have
272    * been converted).
273    *
274    * The final position will be determined utilizing
275    * empirically defined confidence factors for each method.
276    * The final storm position will be returned along with a
277    * position determination flag.
278    *
279    * @param InputLatitudePosition Storm center latitude.
280    * @param InputLongitudePosition Storm center longitude.
281    *
282    * Outputs : Latitude_Return            - final storm center latitude position
283    *           Longitude_Return           - final storm center longitude position
284    *           PositioningMethodID_Return - method used to derive storm location
285    *                      0-error
286    *                      1-interpolation of operational forecast
287    *                      2-Laplacian analysis (not used anymore)
288    *                      3-Warm Spot location
289    *                      4-10^ log spiral analysis
290    *                      5-Combo method of spiral and ring analyses
291    *                      6-linear extrapolation from prior locations
292    * Return  : Error flag = 0
293    */
294   public static double[] AutoMode2(double InputLatitudePosition,
295                                    double InputLongitudePosition)
296       throws IOException
297   {
298      int XInc;
299      int YInc;
300      double SpiralCenterLatitude = -99.9;
301      double SpiralCenterLongitude = -99.9;
302      double SpiralCenterScore = -1.0;
303      double RingFitLatitude = -99.9;
304      double RingFitLongitude = -99.9;
305      double RingFitScore = -1.0;
306      double FinalAutoLatitude = -99.9;
307      double FinalAutoLongitude = -99.9;
308      double FinalAutoScore = -99.9;
309      int FinalAutoFixMethod = 1;
310      boolean DatelineCrossTF = false;
311      boolean DoRemappingTF = true;
312
313      /* NEED TO MULTIPLY LONGITUDES BY -1.0 SINCE ADT ROUTINES WERE DEVELOPED ON McIDAS-X,
314       * which uses West/East Longitudes of +1.0/-1.0.  McV is reversed.
315       */
316      InputLongitudePosition = -1.0 * InputLongitudePosition;  // flip for ADT routines
317      for (YInc = 0; YInc < Data.IRData_NumberRows; YInc++ ) {
318         for (XInc = 0; XInc < Data.IRData_NumberColumns; XInc++ ) {
319             Data.IRData_Longitude[YInc][XInc] = (float)-1.0 * Data.IRData_Longitude[YInc][XInc];
320         }
321      }
322
323      if(DoRemappingTF) {
324         System.out.printf("REMAPPING DATA\n");
325         IRData_Remap_NumberRows = Data.IRData_NumberRows;
326         IRData_Remap_NumberColumns = Data.IRData_NumberColumns;
327
328         double LongitudeIncrement = Data.IRData_Longitude[0][0] - Data.IRData_Longitude[0][1];
329         double LatitudeIncrement = Data.IRData_Latitude[0][0] - Data.IRData_Latitude[1][1];
330
331         if(abs(LongitudeIncrement-LatitudeIncrement)<0.001) {
332            System.out.printf("already remapped\n");
333            /* data is already remapped */
334            /* crosses dateline check */
335            int XSizeMax = Data.IRData_NumberColumns-1;
336            int YSizeMax = Data.IRData_NumberRows-1;
337            double NWCornerLongitude = Data.IRData_Longitude[0][0];
338            double NECornerLongitude = Data.IRData_Longitude[0][XSizeMax];
339            double SWCornerLongitude = Data.IRData_Longitude[YSizeMax][0];
340            double SECornerLongitude = Data.IRData_Longitude[YSizeMax][XSizeMax];
341            /* if((NWCornerLongitude<NECornerLongitude)||(SWCornerLongitude<SECornerLongitude)) {
342               DatelineCrossTF = true;
343            } */
344            if((NWCornerLongitude>NECornerLongitude)||(SWCornerLongitude>SECornerLongitude)) {
345                DatelineCrossTF = true;
346                NWCornerLongitude = NWCornerLongitude+360.0;
347
348            }
349            for(YInc=0; YInc< Data.IRData_NumberRows; YInc++) {
350               for(XInc=0; XInc< Data.IRData_NumberColumns; XInc++) {
351                  IRData_Remap_Longitude[YInc][XInc] = Data.IRData_Longitude[YInc][XInc];
352                  /* check for dateline crossing */
353                  if (DatelineCrossTF&&(IRData_Remap_Longitude[YInc][XInc]<0.0)) {
354                     IRData_Remap_Latitude[YInc][XInc] = Data.IRData_Latitude[YInc][XInc]+360.0;
355                  }
356                  IRData_Remap_Latitude[YInc][XInc] = Data.IRData_Latitude[YInc][XInc];
357                  IRData_Remap_Temperature[YInc][XInc] = Data.IRData_Temperature[YInc][XInc];
358               }
359            }
360            IRData_Remap_NumberColumns = Data.IRData_NumberColumns;
361            IRData_Remap_NumberRows = Data.IRData_NumberRows;
362         } else {
363            /* remap data to rectilinear projection */
364            RemapData( );
365            System.out.printf("COMPLETED REMAPPING DATA\n");
366         }
367    
368         /* perform spiral analysis to determine storm center location */
369         double SpiralReturn[] = SpiralCenterLowRes(InputLatitudePosition,InputLongitudePosition);
370         SpiralCenterLatitude = SpiralReturn[0];
371         SpiralCenterLongitude = SpiralReturn[1];
372         SpiralCenterScore = SpiralReturn[2];
373         
374         /* System.out.printf("Spiral Score : lat=%f  lon=%f  score=%f\n",SpiralCenterLatitude,
375          *                  SpiralCenterLongitude,SpiralCenterScore);
376          */
377            
378         /* redefine first guess for ring analysis as input storm location point */
379         double RingFitFirstGuessLatitude = SpiralCenterLatitude;
380         double RingFitFirstGuessLongitude = SpiralCenterLongitude;
381    
382         /* calculate Moat Mask for false eye check */
383         double MoatMaskTempThreshold = 237.0;
384         double MoatMaskMaxRadiusDegree = 0.50;
385         MoatMaskCalc(MoatMaskTempThreshold,MoatMaskMaxRadiusDegree,1);
386    
387         /* perform ring analysis to determine storm center location */
388         double RingReturn[] = RingFit(RingFitFirstGuessLatitude,RingFitFirstGuessLongitude);
389         RingFitLatitude = RingReturn[0];
390         RingFitLongitude = RingReturn[1];
391         RingFitScore = RingReturn[2];
392         
393         /* System.out.printf("Ring Score   : lat=%f  lon=%f  score=%f \n",
394          *                  RingFitLatitude,RingFitLongitude,RingFitScore);
395          */
396         /*
397          * System.out.printf("fglat=%f fglon=%f  SpiralCenterLatitude=%f SpiralCenterLongitude=%f rglat=%f rglon=%f\n",
398          *                   InputLatitudePosition,InputLongitudePosition,
399          *                   SpiralCenterLatitude,SpiralCenterLongitude,
400          *                   RingFitLatitude,RingFitLongitude);
401          */
402         
403         /* caluculate confidence factor for combined spiral/ring analyses */
404         double ScoresReturn[] = CalcScores(InputLatitudePosition,InputLongitudePosition,
405                                            SpiralCenterLatitude,SpiralCenterLongitude,SpiralCenterScore,
406                                            RingFitLatitude,RingFitLongitude,RingFitScore);
407         FinalAutoLatitude = ScoresReturn[0];
408         FinalAutoLongitude = ScoresReturn[1];
409         FinalAutoScore = ScoresReturn[2];
410         FinalAutoFixMethod = (int)ScoresReturn[3];
411         /*
412          * System.out.printf("SPIRAL CENTER : lat=%f  lon=%f  SCORE=%f  \n",
413          *                    SpiralCenterLatitude,SpiralCenterLongitude,SpiralCenterScore);
414          * System.out.printf("RING   CENTER : lat=%f  lon=%f  SCORE=%f  \n",
415          *                    RingFitLatitude,RingFitLongitude,RingFitScore);
416          * System.out.printf("FINAL  CENTER : lat=%f  lon=%f  SCORE=%f  \n",
417          *                    FinalAutoLatitude,FinalAutoLongitude,FinalAutoScore);
418          */
419      } else {
420         SpiralCenterLatitude = -99.9;
421         SpiralCenterLongitude = -99.9;
422         SpiralCenterScore = -1.0;
423         RingFitLatitude = -99.9;
424         RingFitLatitude = -99.9;
425         RingFitScore = -1.0;
426         FinalAutoLatitude = -99.9;
427         FinalAutoLongitude = -99.9;
428         FinalAutoScore = -99.9;
429         FinalAutoFixMethod = 1;
430      }  
431
432      int PositioningMethodID = History.IRCurrentRecord.autopos;
433    
434      double LocationReturn[] = PickFinalLocation(PositioningMethodID,
435                                                  InputLatitudePosition,
436                                                  InputLongitudePosition,
437                                                  FinalAutoLatitude,
438                                                  FinalAutoLongitude,
439                                                  FinalAutoScore,
440                                                  FinalAutoFixMethod);
441      
442      double FinalLatitude = LocationReturn[0];
443      double FinalLongitude = LocationReturn[1];
444      PositioningMethodID = (int)LocationReturn[3];
445    
446      if (FinalLongitude > 180.0) {
447         FinalLongitude = FinalLongitude-360.0;
448      }
449    
450      /* System.out.printf("TAC CENTER    :  lat=%f  lon=%f \n",
451                         InputLatitudePosition,InputLongitudePosition); */
452      /* System.out.printf("AUTO CENTER   :  lat=%f  lon=%f  SCORE=%f  METHOD=%d\n",
453                         FinalAutoLatitude,FinalAutoLongitude,
454                         FinalAutoScore,FinalAutoFixMethod); */
455      /* System.out.printf("FINAL LOCATION:  lat=%f  lon=%f  SCORE=%f  METHOD=%d\n",
456                         FinalLatitude,FinalLongitude,
457                         FinalScoreValue,PositioningMethodID); */
458      
459      /* need to flip final Longitude value back to McV format by multiplying by -1.0 */
460      FinalLongitude = -1.0 * FinalLongitude;
461
462      return new double[] { FinalLatitude, FinalLongitude, (double)PositioningMethodID };
463   }
464
465   private static double[] SpiralCenterLowRes(double InputLatitude, double InputLongitude) {
466
467      double[][] IRData_NormCoord_Latitude = new double[200][200];
468      double[][] IRData_NormCoord_Longitude = new double[200][200];
469      double[][] IRData_NormCoord_Temperature = new double[200][200];
470
471      int XInc,YInc;
472      int IRData_NormCoord_NumberRows;
473      int IRData_NormCoord_NumberColumns;
474      int NumPoints;
475      int RingWidthLocal = (int)RING_WIDTH;
476
477      double XOff,YOff;
478      double SearchRadius = 0.0;
479      double CrossScoreSum = 0.0;
480      double ProxyValueX = 0.0;
481      double ProxyValueY = 0.0;
482      double DenominatorValue = 0.0;
483      double SpiralValueX = 0.0;
484      double SpiralValueY = 0.0;
485      double RawCrossScore = 0.0;
486      double CrossScoreClean = 0.0;
487      double SpiralMeanCrossScore = 0.0;
488      double SpiralMeanCrossMaxScore = 0.0;
489      double SpiralMeanCrossMaxYLocation = 0.0;
490      double SpiralMeanCrossMaxXLocation = 0.0;
491
492      double SpiralCenterLatitude = -99.99;
493      double SpiralCenterLongitude = -999.99;
494      double SpiralCenterScore = 0.0;
495
496      double Alpha = 5.0*PI/180.0;
497      double OutsideFactor = 0.62;
498
499      double OuterSearchRadiusDegree = 1.75;
500      double CourseGridSpacingDegree = 0.2;
501      double FineGridSpacingDegree = 0.1;
502      double FilterDiscSize = pow(OuterSearchRadiusDegree+(2.0*FineGridSpacingDegree),2);
503      double AlphaPOWP1 = 1.0+pow(Alpha,2);
504
505      int ImageResolution = (int) Data.GetCurrentImageResolution();
506      int IncAddVal = (ImageResolution>RingWidthLocal) ? 1 : (RingWidthLocal-ImageResolution+1);
507
508      double SignFactor = abs(InputLatitude)/InputLatitude;
509
510      YInc = (IRData_Remap_NumberRows)-1;
511      if((InputLongitude<0.0)&&((IRData_Remap_Longitude[0][0]>180.0)||
512         (IRData_Remap_Longitude[YInc][0]>180.0))) {
513         /*
514          * System.out.printf("DATELINE CROSS... changing InputLongitude from %f",InputLongitude);
515          */
516         InputLongitude = InputLongitude+360.0;
517      }
518
519      for(YInc=0;YInc<IRData_Remap_NumberRows;YInc++) {
520         for(XInc=0;XInc<IRData_Remap_NumberColumns;XInc++) {
521            /* compute normalized coordinate system arrays */
522            IRData_NormCoord_Longitude[YInc][XInc] = (IRData_Remap_Longitude[YInc][XInc]-InputLongitude)*
523                                                 (cos(PI*InputLatitude/180.0));
524            IRData_NormCoord_Latitude[YInc][XInc] = IRData_Remap_Latitude[YInc][XInc]-InputLatitude;
525            IRData_NormCoord_Temperature[YInc][XInc] = IRData_Remap_Temperature[YInc][XInc];
526            /*
527             * System.out.printf("YInc=%d XInc=%d  lat=%f lon=%f  lat=%f lon=%f  temp=%f\n",
528             *                    YInc,XInc,
529             *                    IRData_Remap_Latitude[YInc][XInc],
530             *                    IRData_Remap_Longitude[YInc][XInc],
531             *                    IRData_NormCoord_Latitude[YInc][XInc],
532             *                    IRData_NormCoord_Longitude[YInc][XInc],
533             *                    IRData_Remap_Temperature[YInc][XInc]);
534             */
535         }
536      }
537      IRData_NormCoord_NumberColumns = IRData_Remap_NumberColumns;
538      IRData_NormCoord_NumberRows = IRData_Remap_NumberRows;
539
540      /* determine lat/lon grid increment */
541      /* W to E gradient */
542      double LongitudeIncrement = abs(IRData_NormCoord_Longitude[0][0]-
543                                           IRData_NormCoord_Longitude[0][1]);
544      /* N to S gradient */
545      double LatitudeIncrement = abs(IRData_NormCoord_Latitude[0][0]-
546                                          IRData_NormCoord_Latitude[1][0]);
547      /*
548       * This is to determine longitude multiplier factor... original routines
549       * were developed using negative, but McIDAS is positive in WH.  So if
550       * LongitudeMultFactor is negative, we know (assume) we are using non-McIDAS
551       * input imagery/values, otherwise make LongitudeMultFactor positive.
552       * This all assumes that image is loaded from NW to SE
553       */
554      double LongitudeMultFactor = IRData_NormCoord_Longitude[0][0]-IRData_NormCoord_Longitude[0][1];
555      LongitudeMultFactor = (LongitudeMultFactor<0.0) ? 1.0 : -1.0;
556      /*
557       * System.out.printf("LatitudeIncrement=%f  LongitudeIncrement=%f LongitudeMultFactor=%f\n",
558       *                    LatitudeIncrement, LongitudeIncrement,LongitudeMultFactor);
559       */
560      
561    
562      /* calculate gradient field */
563      Gradient(IRData_NormCoord_Temperature,IRData_NormCoord_NumberColumns,IRData_NormCoord_NumberRows,
564               LongitudeIncrement,LatitudeIncrement);
565
566      for(YInc=0;YInc<IRData_NormCoord_NumberRows-1;YInc++) {
567         for(XInc=0;XInc<IRData_NormCoord_NumberColumns-1;XInc++) {
568            /*
569             * System.out.printf("YInc=%d XInc=%d  remap:lat=%f lon=%f  normcoord:lat=%f lon=%f NSTempGradientArray=%f EWTempGrad=%f\n",
570             *                 YInc,XInc,
571             *                 IRData_Remap_Latitude[YInc][XInc],IRData_Remap_Longitude[YInc][XInc],
572             *                 IRData_NormCoord_Latitude[YInc][XInc],IRData_NormCoord_Longitude[YInc][XInc],
573             *                 NSTempGradientArray[YInc][XInc],EWTempGradientArray[YInc][XInc]);
574             */
575   
576            double GradientOriginalMag = sqrt(pow(NSTempGradientArray[YInc][XInc],2)+
577                                              pow(EWTempGradientArray[YInc][XInc],2));
578            double GradientLOGMag = log(1.0+GradientOriginalMag);
579            double GradientLOGReduction = 0.0;
580            if(GradientLOGMag!=0.0) {
581               GradientLOGReduction = GradientLOGMag/GradientOriginalMag;
582            }
583            NSTempGradientArray[YInc][XInc] = GradientLOGReduction*NSTempGradientArray[YInc][XInc];
584            EWTempGradientArray[YInc][XInc] = GradientLOGReduction*EWTempGradientArray[YInc][XInc];
585         }
586      }
587    
588      /* COURSE GRID */
589      /* calculate cross product score at each grid point */
590      /* XInc/YInc are "starting point" coordinates */
591      SpiralMeanCrossMaxScore = -99.0;
592      double SearchRadiusMaximum = pow(OuterSearchRadiusDegree+((2.0*CourseGridSpacingDegree)/3.0),2);
593      for(XOff=-OuterSearchRadiusDegree;XOff<=OuterSearchRadiusDegree;XOff=XOff+CourseGridSpacingDegree) {
594         for(YOff=-OuterSearchRadiusDegree;YOff<=OuterSearchRadiusDegree;YOff=YOff+CourseGridSpacingDegree) {
595            /* XOff/YOff are offset coordinates from "starting point" */
596            SearchRadius = pow(XOff,2)+pow(YOff,2);
597            if(SearchRadius<=SearchRadiusMaximum) {
598               CrossScoreSum=0.0;
599               NumPoints=0;
600               for(YInc=1;YInc<IRData_NormCoord_NumberRows-2;YInc=YInc+IncAddVal) {
601                  for(XInc=1;XInc<IRData_NormCoord_NumberColumns-2;XInc=XInc+IncAddVal) {
602                     SearchRadius = pow(IRData_NormCoord_Longitude[YInc][XInc],2)+
603                                    pow(IRData_NormCoord_Latitude[YInc][XInc],2);
604                     if(SearchRadius<FilterDiscSize) {
605                        ProxyValueX = LongitudeMultFactor*IRData_NormCoord_Longitude[YInc][XInc]-XOff;
606                        ProxyValueY = IRData_NormCoord_Latitude[YInc][XInc]-YOff;
607                        DenominatorValue = sqrt((AlphaPOWP1)*(pow(ProxyValueX,2)+pow(ProxyValueY,2)));
608                        SpiralValueX = ((Alpha*ProxyValueX)+(SignFactor*ProxyValueY))/DenominatorValue;
609                        SpiralValueY = ((Alpha*ProxyValueY)-(SignFactor*ProxyValueX))/DenominatorValue;
610                        RawCrossScore = (SpiralValueX*NSTempGradientArray[YInc][XInc])-
611                                        (SpiralValueY* EWTempGradientArray[YInc][XInc]);
612                        CrossScoreClean = max(0.0,-RawCrossScore)+(OutsideFactor*max(0.0,RawCrossScore));
613                        CrossScoreSum = CrossScoreSum+CrossScoreClean;
614                        /*
615                         * System.out.printf("%d %d : lat=%f lon=%f  xoff=%f yoff=%f ",
616                         *                   "ProxyValueX=%f ProxyValueY=%f ",
617                         *                   "SpiralValueX=%f SpiralValueY=%f ",
618                         *                   "rawScore=%f clean=%f\n",
619                         *                    IRData_Remap_Latitude[YInc][XInc],
620                         *                    IRData_Remap_Longitude[YInc][XInc],
621                         *                    XOff,YOff,ProxyValueX,ProxyValueY,
622                         *                    SpiralValueX,SpiralValueY,
623                         *                    RawCrossScore,CrossScoreClean);
624                         */
625                        NumPoints++;
626                     }
627                  }
628               }
629               /* calculate mean of all values in CrossScore array */
630               SpiralMeanCrossScore = CrossScoreSum/(double)NumPoints;
631               /* store location of maximum score position */
632               if(SpiralMeanCrossScore>SpiralMeanCrossMaxScore) {
633                  SpiralMeanCrossMaxScore = SpiralMeanCrossScore;
634                  SpiralMeanCrossMaxYLocation = YOff;
635                  SpiralMeanCrossMaxXLocation = XOff;
636               }
637            }
638         }
639      }
640     
641       /* System.out.printf("course grid : y=%f x=%f max=%f\n",
642                          SpiralMeanCrossMaxYLocation,SpiralMeanCrossMaxXLocation,
643                          SpiralMeanCrossMaxScore); */
644      
645    
646      /* FINE GRID */
647      int SpiralCount = 1;
648      double FineGridXMinimum = SpiralMeanCrossMaxXLocation-CourseGridSpacingDegree;
649      double FineGridXMaximum = SpiralMeanCrossMaxXLocation+CourseGridSpacingDegree;
650      double FineGridYMimimum = SpiralMeanCrossMaxYLocation-CourseGridSpacingDegree;
651      double FineGridYMaximum = SpiralMeanCrossMaxYLocation+CourseGridSpacingDegree;
652      SpiralMeanCrossMaxScore=-99.0;
653      for(XOff=FineGridXMinimum;XOff<=FineGridXMaximum;XOff=XOff+FineGridSpacingDegree) {
654         for(YOff=FineGridYMimimum;YOff<=FineGridYMaximum;YOff=YOff+FineGridSpacingDegree) {
655            /* XOff/YOff are offset coordinates from "starting point" */
656            CrossScoreSum = 0.0;
657            NumPoints = 0;
658            for(YInc=1;YInc<IRData_NormCoord_NumberRows-2;YInc++) {
659               for(XInc=1;XInc<IRData_NormCoord_NumberColumns-2;XInc++) {
660                  SearchRadius = Math.pow(IRData_NormCoord_Longitude[YInc][XInc],2)+
661                                 Math.pow(IRData_NormCoord_Latitude[YInc][XInc],2);
662                  if(SearchRadius<FilterDiscSize) {
663                      ProxyValueX = LongitudeMultFactor*IRData_NormCoord_Longitude[YInc][XInc]-XOff;
664                      ProxyValueY = IRData_NormCoord_Latitude[YInc][XInc]-YOff;
665                      DenominatorValue = Math.sqrt((AlphaPOWP1)*(Math.pow(ProxyValueX,2)+Math.pow(ProxyValueY,2)));
666                      SpiralValueX = (Alpha*ProxyValueX+(SignFactor*ProxyValueY))/DenominatorValue;
667                      SpiralValueY = (Alpha*ProxyValueY-(SignFactor*ProxyValueX))/DenominatorValue;
668                      RawCrossScore = (SpiralValueX*NSTempGradientArray[YInc][XInc])-
669                                      (SpiralValueY*EWTempGradientArray[YInc][XInc]);
670                      CrossScoreClean = Math.max(0.0,-RawCrossScore)+(OutsideFactor*Math.max(0.0,RawCrossScore));
671                      CrossScoreSum = CrossScoreSum+CrossScoreClean;
672                      NumPoints++;
673                   }
674                }
675             }
676             /* calculate mean of all values in CrossScore array */
677             SpiralMeanCrossScore = CrossScoreSum/(double)NumPoints;
678             /* store location of maximum score position */
679             if(SpiralMeanCrossScore>SpiralMeanCrossMaxScore) {
680                SpiralMeanCrossMaxScore = SpiralMeanCrossScore;
681                SpiralMeanCrossMaxYLocation = YOff;
682                SpiralMeanCrossMaxXLocation = XOff;
683             }
684             SpiralCenterAnalysisField[SpiralCount][0] = SpiralMeanCrossScore;
685             SpiralCenterAnalysisField[SpiralCount][1] = (double)YOff;
686             SpiralCenterAnalysisField[SpiralCount][2] = (double)XOff;
687             SpiralCenterAnalysisField[0][0] = (double)SpiralCount;
688             SpiralCount++;
689         }
690      }
691      
692      /* System.out.printf("fine grid : y=%f x=%f max=%f\n",
693                          SpiralMeanCrossMaxYLocation,SpiralMeanCrossMaxXLocation,SpiralMeanCrossMaxScore); */
694      
695    
696      SpiralCenterAnalysisField[0][1] = 0.0;
697      SpiralCenterAnalysisField[0][2] = 0.0;
698    
699      /* determine lat/lon point from x/y coordinates */
700      SpiralCenterLatitude = SpiralMeanCrossMaxYLocation+InputLatitude;
701      SpiralCenterLongitude = ((LongitudeMultFactor*SpiralMeanCrossMaxXLocation)/
702                                     (Math.cos(PI*InputLatitude/180.0)))+InputLongitude;
703      SpiralCenterScore = SpiralMeanCrossMaxScore;
704      for(YInc=1;YInc<=(int)SpiralCenterAnalysisField[0][0];YInc++) {
705         SpiralCenterAnalysisField[YInc][1] = SpiralCenterAnalysisField[YInc][1]+InputLatitude;
706         SpiralCenterAnalysisField[YInc][2] = ((LongitudeMultFactor*SpiralCenterAnalysisField[YInc][2])/
707                                               (Math.cos(PI*InputLatitude/180.0)))+InputLongitude;
708      }
709
710      return new double[] { SpiralCenterLatitude, SpiralCenterLongitude, SpiralCenterScore };
711   }
712
713   private static void Gradient(double TemperatureInputArray[][], int ElementXNumber, int LineYNumber,
714                                double LongitudeIncrement, double LatitudeIncrement) {
715
716      int XInc,YInc;
717
718      /* initialize arrays */
719      for(YInc=0;YInc<LineYNumber-1;YInc++) {
720         for(XInc=0;XInc<ElementXNumber-1;XInc++) {
721            NSTempGradientArray[YInc][XInc] = 0.0;
722            EWTempGradientArray[YInc][XInc] = 0.0;
723         }
724      }
725
726      for(YInc=1;YInc<LineYNumber-1;YInc++) {
727         for(XInc=1;XInc<ElementXNumber-1;XInc++) {
728            /* determine N-S gradient at point */
729            NSTempGradientArray[YInc][XInc] = (TemperatureInputArray[YInc-1][XInc]-
730                                               TemperatureInputArray[YInc+1][XInc])/
731                                               (2.0*LatitudeIncrement);
732            /* determine E-W gradient at point */
733            EWTempGradientArray[YInc][XInc] = (TemperatureInputArray[YInc][XInc+1]-
734                                               TemperatureInputArray[YInc][XInc-1])/
735                                               (2.0*LongitudeIncrement);
736         }
737      }
738
739   }
740
741   private static double[] RingFit(double RingFitFirstGuessLatitude, double RingFitFirstGuessLongitude) {
742
743      int XInc,YInc,ZInc;
744      int CircleFilterRadii;
745      int CircleFilterPoints;
746
747      int RingWidthLocal = (int)RING_WIDTH;
748      int XLocation = 0;
749      int YLocation = 0;
750      int CircleFilterXLocation = 0;
751      int CircleFilterYLocation = 0;
752
753      double[][] GradientFieldXDirection = new double[200][200];
754      double[][] GradientFieldYDirection = new double[200][200];
755
756      double RingFitSearchRadiusDegree = 0.75;
757      double RingFitMinRadiusDegree = 0.06;
758      double RingFitMaxRadiusDegree = 0.40;
759      double CircleFilterDistance = 0.0;
760
761      double DotScoreMaximum = -99999.0;
762      double DotProductXDirection = -999.9;
763      double DotProductYDirection = -999.9;
764      double DotProductXYValue = -999.9;
765      double SignFactor = 0.0;
766      double DotScoreValue = 0.0;
767      double DotScoreValueSum = 0.0;
768      double DotScoreFinal = 0.0;
769      double NANAdjustmentValue = 0.0;
770
771      int ImageResolution = (int) Data.GetCurrentImageResolution();
772      int IncAddVal = (ImageResolution>RingWidthLocal) ? 1 : (RingWidthLocal-ImageResolution+1);
773
774      /* derive values */
775      double DegreesPerPixel = Math.abs(IRData_Remap_Latitude[0][0]-IRData_Remap_Latitude[1][0]);
776      int RingFitSearchRadiusPixel = (int)Math.round(RingFitSearchRadiusDegree/DegreesPerPixel);
777      int RingFitMinRadiusPix = (int)Math.max(2.0,Math.round(RingFitMinRadiusDegree/DegreesPerPixel));
778      int MaximumRadiusSizePixels = (int)Math.round(RingFitMaxRadiusDegree/DegreesPerPixel);
779
780      double LongitudeIncrement = 1.0;
781      double LatitudeIncrement = -1.0;
782
783      
784      /* System.out.printf("RingFit: %d %f %d %d %d\n",IncAddVal,DegreesPerPixel,
785                                 RingFitSearchRadiusPixel,RingFitMinRadiusPix,
786                                 MaximumRadiusSizePixels); */
787      /*
788       * for(YInc=0;YInc<IRData_Remap_NumberRows-1;YInc++) {
789       *    for(XInc=0;XInc<IRData_Remap_NumberColumns-1;XInc++) {
790       *       System.out.printf("Y=%d X=%d lat=%f lon=%f\n",YInc,XInc,IRData_Remap_Latitude[YInc][XInc],IRData_Remap_Longitude[YInc][XInc]);
791       *    }
792       *  }
793       */
794
795      /* NEED TO PASS BACK BOTH 2D ARRAYS FROM GRADIENT AND CIRCLEFILT */
796      /* calculate gradient field */
797      Gradient(IRData_Remap_Temperature,IRData_Remap_NumberColumns,IRData_Remap_NumberRows,
798               LongitudeIncrement,LatitudeIncrement);
799
800      for(YInc=0;YInc<IRData_Remap_NumberRows-1;YInc++) {
801         for(XInc=0;XInc<IRData_Remap_NumberColumns-1;XInc++) {
802            GradientFieldYDirection[YInc][XInc] = NSTempGradientArray[YInc][XInc];
803            GradientFieldXDirection[YInc][XInc] = EWTempGradientArray[YInc][XInc];
804            /* System.out.printf("GRADIENT  Y=%d X=%d Xdir=%f Ydir=%f  ",YInc,XInc,GradientFieldXDirection[YInc][XInc],GradientFieldYDirection[YInc][XInc]); */
805            /* System.out.printf(" lat=%f lon=%f\n",IRData_Remap_Latitude[YInc][XInc],IRData_Remap_Longitude[YInc][XInc]); */
806         }
807      }
808      /*
809       * System.out.printf("RingFit : RingFitFirstGuessLatitude=%f RingFitFirstGuessLongitude=%f\n",
810       *                    RingFitFirstGuessLatitude,RingFitFirstGuessLongitude);
811       */
812
813      /* make matricies of row and column numbers */
814      int Lalo2IndsReturn[] = Lalo2IndsFloat(RingFitFirstGuessLatitude,RingFitFirstGuessLongitude,
815                                             IRData_Remap_Latitude,IRData_Remap_Longitude,
816                                             IRData_Remap_NumberColumns,IRData_Remap_NumberRows);
817
818      int FirstGuessXLocation = Lalo2IndsReturn[0];
819      int FirstGuessYLocation = Lalo2IndsReturn[1];
820
821      /*
822       * System.out.printf("RingFit : FirstGuessXLocation=%d FirstGuessYLocation=%d\n",
823       *                    FirstGuessXLocation,FirstGuessYLocation);
824       */
825    
826      /* initialize circle/ring filter arrays */
827      /* radius in pixels */
828
829      for(CircleFilterRadii=0;CircleFilterRadii<100;CircleFilterRadii++) {
830         /* number of points on circle at radius */
831         for(CircleFilterPoints=0;CircleFilterPoints<2000;CircleFilterPoints++) {
832            CircleFilterRowArray[CircleFilterRadii][CircleFilterPoints] = 0;
833            CircleFilterColumnArray[CircleFilterRadii][CircleFilterPoints] = 0;
834         }
835      }
836
837      /*
838       * determine digital pixel coordinates for ring analysis for
839       * different radii sizes
840       */
841      for(CircleFilterRadii=RingFitMinRadiusPix;CircleFilterRadii<=MaximumRadiusSizePixels;CircleFilterRadii++) {
842         CircleFilt(CircleFilterRadii);
843      }
844    
845      /* search image box */
846      ZInc = 1;
847      /* develop the accumulator */
848      for(CircleFilterRadii=RingFitMinRadiusPix;CircleFilterRadii<MaximumRadiusSizePixels;CircleFilterRadii++) {
849         int CircleFilterPointCount = CircleFilterRowArray[CircleFilterRadii][0];
850         /* determine each main point in analysis disc */
851         for(XInc=1;XInc<IRData_Remap_NumberColumns-1;XInc=XInc+IncAddVal) {
852            for(YInc=1;YInc<IRData_Remap_NumberRows-1;YInc=YInc+IncAddVal) {
853               /*
854                * System.out.printf("XInc=%d YInc=%d  %d %d: %d\n",XInc,YInc,
855                *                    FirstGuessXLocation,FirstGuessYLocation,RingFitSearchRadiusPixel);
856                */
857
858               CircleFilterDistance = (double)((XInc-FirstGuessXLocation)*(XInc-FirstGuessXLocation))+
859                                              ((YInc-FirstGuessYLocation)*(YInc-FirstGuessYLocation));
860               /*
861                * System.out.printf("XInc=%d YInc=%d CircleFilterDistance=%f RingFitSearchRadiusPixel=%d\n",
862                *                    XInc,YInc,CircleFilterDistance,RingFitSearchRadiusPixel);
863                */
864               if(CircleFilterDistance<=(double)(RingFitSearchRadiusPixel*RingFitSearchRadiusPixel)) {
865                  /*
866                   * if main point (YInc,XInc) is in disc, calculate dotproduct
867                   * for each subpoint on ring around main point
868                   */
869                  DotScoreValueSum = 0.0;
870                  int NANCount = 0;
871                  boolean FoundMoatRegionTF = false;
872                  /*
873                   * System.out.printf("XInc=%d YInc=%d  CircleFilterPointCount=%d %d\n",
874                   *                    XInc,YInc,CircleFilterPointCount,CircleFilterRadii);
875                   */
876                  for(CircleFilterPoints=1;CircleFilterPoints<=CircleFilterPointCount;CircleFilterPoints++) {
877                     /*
878                      * System.out.printf("CircleFilterPoints=%d CircleFilterRadii=%d XInc=%d YInc=%d\n",
879                      *                    CircleFilterPoints,CircleFilterRadii,XInc,YInc);
880                      */
881                     CircleFilterXLocation = XInc+CircleFilterColumnArray[CircleFilterRadii][CircleFilterPoints];
882                     CircleFilterYLocation = YInc+CircleFilterRowArray[CircleFilterRadii][CircleFilterPoints];
883                     /*
884                      * System.out.printf("%d : %d %d --  %d : %d %d ",
885                      *                    CircleFilterYLocation,YInc,
886                      *                    CircleFilterRowArray[CircleFilterRadii][CircleFilterPoints],
887                      *                    CircleFilterXLocation,XInc,
888                      *                    CircleFilterColumnArray[CircleFilterRadii][CircleFilterPoints]);
889                      */
890                     if((CircleFilterXLocation<1)||(CircleFilterYLocation<1)) {
891                        DotProductXDirection = -999.9;
892                        DotProductYDirection = -999.9;
893                     } else {
894                        /*
895                         * System.out.printf("I=%d x=%d  J=%d y=%d  : %b %b\n",
896                         *                    CircleFilterXLocation,IRData_Remap_NumberColumns,
897                         *                    CircleFilterYLocation,IRData_Remap_NumberRows,
898                         *                    (CircleFilterXLocation>=IRData_Remap_NumberColumns),
899                         *                    (CircleFilterYLocation>=IRData_Remap_NumberRows));
900                         */
901                        if((CircleFilterXLocation>=(IRData_Remap_NumberColumns-1))||
902                           (CircleFilterYLocation>=(IRData_Remap_NumberRows-1))) {
903                           DotProductXDirection = -999.9;
904                           DotProductYDirection = -999.9;
905                        } else {
906                           if(MoatMaskFlagField[CircleFilterYLocation][CircleFilterXLocation]==1) {
907                              FoundMoatRegionTF = true;
908                           }
909                           /* System.out.printf(" mask=%d  FoundMoatRegionTF=%b\n",MoatMaskFlagField[CircleFilterYLocation][CircleFilterXLocation],FoundMoatRegionTF); */
910                           if(FoundMoatRegionTF) {
911                              DotProductXDirection = -999.9;
912                              DotProductYDirection = -999.9;
913                           } else {
914                              DotProductXDirection=
915                              ((double)CircleFilterRowArray[CircleFilterRadii][CircleFilterPoints]/
916                               (double)CircleFilterRadii)*
917                               GradientFieldYDirection[CircleFilterYLocation][CircleFilterXLocation];
918                               /*
919                                * System.out.printf(" XX | %d %d j=%d i=%d j=%d i=%d  GradientFieldYDirection=%f  GradientFieldXDirection=%f | \n",
920                                *                    CircleFilterPoints,CircleFilterRadii,
921                                *                    CircleFilterRowArray[CircleFilterRadii][CircleFilterPoints],
922                                *                    CircleFilterColumnArray[CircleFilterRadii][CircleFilterPoints],
923                                *                    CircleFilterYLocation,
924                                *                    CircleFilterXLocation,
925                                *                    GradientFieldYDirection[CircleFilterYLocation][CircleFilterXLocation],
926                                *                    GradientFieldXDirection[CircleFilterYLocation][CircleFilterXLocation]);
927                                */
928                              DotProductYDirection=
929                               ((double)CircleFilterColumnArray[CircleFilterRadii][CircleFilterPoints]/
930                                (double)CircleFilterRadii)*
931                                 GradientFieldXDirection[CircleFilterYLocation][CircleFilterXLocation];
932                           }
933                        }
934                     }
935                     if((DotProductXDirection<-999.0)||(DotProductYDirection<-999.0)) {
936                        NANCount++;
937                     } else {
938                        DotProductXYValue = DotProductXDirection+DotProductYDirection;
939                        if(DotProductXYValue==0.0) {
940                           SignFactor = 0.0;
941                        } else {
942                           /* return -1/+1 for -/+ value */
943                           SignFactor = Math.abs(DotProductXYValue)/DotProductXYValue;
944                        }
945                        DotScoreValue = SignFactor*(Math.log(1.0+Math.abs(DotProductXYValue)));
946                        DotScoreValueSum = DotScoreValueSum+DotScoreValue;
947                     }
948                     /*
949                      * System.out.printf("dot product X=%f Y=%f XY=%f value=%f sum=%f\n",DotProductXDirection,
950                      *            DotProductYDirection,DotProductXYValue,DotScoreValue,DotScoreValueSum);
951                      */
952                  } /* if indisk */
953                  /* System.out.printf("dot product Final Sum=%f \n",DotScoreValueSum); */
954                  /* check for missing data and adjust DotScoreFinal accordingly */
955                  /* System.out.printf("Y=%d X=%d foundmoat=%b  nancount=%f  xyz=%f\n",YInc,XInc,FoundMoatRegionTF,(double)NANCount,0.575*(double)CircleFilterPointCount); */
956                  if(FoundMoatRegionTF||((double)NANCount)>(0.575*(double)CircleFilterPointCount)) {
957                     DotScoreFinal = 0.0;
958                  } else {
959                     NANAdjustmentValue = (double)CircleFilterPointCount/(double)(CircleFilterPointCount-NANCount);
960                     DotScoreFinal = -NANAdjustmentValue*DotScoreValueSum/Math.sqrt((double)CircleFilterPointCount);
961                  }
962                  /* System.out.printf("XX final=%f maximum=%f ",DotScoreFinal,DotScoreMaximum); */
963                  /* System.out.printf("    circlefilterpointcount=%d sqrt=%f ",CircleFilterPointCount,Math.sqrt((double)CircleFilterPointCount)); */
964                  if(DotScoreFinal>DotScoreMaximum) {
965                     DotScoreMaximum = DotScoreFinal;
966                     XLocation = XInc;
967                     YLocation = YInc;
968                     /* System.out.printf("   xloc=%d yloc=%d MaxRadSize=%d",XInc,YInc,MaximumRadiusSize); */
969                  }
970                  /* System.out.printf("\n"); */
971                  RingScoreAnalysisField[ZInc][0] = DotScoreFinal;
972                  RingScoreAnalysisField[ZInc][1] = IRData_Remap_Latitude[YInc][XInc];
973                  RingScoreAnalysisField[ZInc][2] = IRData_Remap_Longitude[YInc][XInc];
974                  RingScoreAnalysisField[0][0] = (double)ZInc;
975                  /*
976                   * System.out.printf("ZZZZ  0=%f 1=%f 2=%f x=%f\n",RingScoreAnalysisField[ZInc][0],RingScoreAnalysisField[ZInc][1],
977                   *                                                 RingScoreAnalysisField[ZInc][2],RingScoreAnalysisField[0][0]);
978                   */
979                  ZInc++;
980               }
981            } /* XInc */
982         } /* YInc */
983      } /* CircleFilterRadii */
984    
985      RingScoreAnalysisField[0][1] = 0.0;
986      RingScoreAnalysisField[0][2] = 0.0;
987    
988      /* make matricies of row and column numbers */
989      /* System.out.printf("inds2lalo: x=%d y=%d \n",XLocation,YLocation); */
990      double Inds2LaloReturn[] = Inds2LaloFloat(XLocation,YLocation,
991                                                IRData_Remap_Latitude,IRData_Remap_Longitude,
992                                                IRData_Remap_NumberColumns,IRData_Remap_NumberRows);
993      double MaxScoreLatitude = Inds2LaloReturn[0];
994      double MaxScoreLongitude = Inds2LaloReturn[1];
995
996      double RingFitLatitude = MaxScoreLatitude;
997      double RingFitLongitude = MaxScoreLongitude;
998      double RingFitMaxScore = DotScoreMaximum;
999
1000      return new double[] { RingFitLatitude, RingFitLongitude, RingFitMaxScore };
1001   }
1002
1003   private static int[] Lalo2IndsFloat( double LatitudeInput, double LongitudeInput,
1004                                        double LatitudeArrayInput[][], double LongitudeArrayInput[][],
1005                                        int ElementXNumber, int LineYNumber) {
1006
1007      double LatitudeMinimum = LatitudeArrayInput[LineYNumber-1][0];
1008      double LatitudeMaximum = LatitudeArrayInput[0][0];
1009      double LongitudeMinimum = LongitudeArrayInput[0][ElementXNumber-1];
1010      double LongitudeMaximum = LongitudeArrayInput[0][0];
1011    
1012      int YAxisPositionReturn = (int)((((double)LineYNumber-1.0)/(LatitudeMaximum-LatitudeMinimum))*
1013                                        (LatitudeMaximum-LatitudeInput));
1014      int XAxisPositionReturn = (int)((((double)ElementXNumber-1.0)/(LongitudeMaximum-LongitudeMinimum))*
1015                                        (LongitudeMaximum-LongitudeInput));
1016
1017      return new int[] { XAxisPositionReturn, YAxisPositionReturn };
1018
1019   } 
1020
1021   private static double[] Inds2LaloFloat(int XAxisPosition, int YAxisPosition,
1022                                          double LatitudeArrayInput[][], double LongitudeArrayInput[][],
1023                                          int ElementXNumber, int LineYNumber) {
1024
1025      double LatitudeMinimum = LatitudeArrayInput[LineYNumber-1][0];
1026      double LatitudeMaximum = LatitudeArrayInput[0][0];
1027      double LongitudeMinimum = LongitudeArrayInput[0][ElementXNumber-1];
1028      double LongitudeMaximum = LongitudeArrayInput[0][0];
1029    
1030      double LongitudeReturn = LongitudeMaximum-((((double)XAxisPosition)/((double)ElementXNumber-1.0))*
1031                                                   (LongitudeMaximum-LongitudeMinimum));
1032      double LatitudeReturn = LatitudeMaximum-((((double)YAxisPosition)/((double)LineYNumber-1.0))*
1033                                                 (LatitudeMaximum-LatitudeMinimum));
1034
1035      return new double [] { LatitudeReturn, LongitudeReturn };
1036
1037   }
1038
1039   private static void CircleFilt(int RingRadiusInput) {
1040
1041      int XInc,YInc;
1042      int PointCount = 1;
1043      int RadiusPlus1 = RingRadiusInput+1;
1044
1045      double ThresholdDifference=0.5*(double)(((RingRadiusInput+1)*(RingRadiusInput+1))-
1046                                               (RingRadiusInput*RingRadiusInput));
1047      for(XInc=-RadiusPlus1;XInc<=RadiusPlus1;XInc++) {
1048         for(YInc=-RadiusPlus1;YInc<=RadiusPlus1;YInc++) {
1049            double DifferenceValue = (double)((YInc*YInc)+(XInc*XInc)-(RingRadiusInput*RingRadiusInput));
1050            if(Math.abs(DifferenceValue)<=ThresholdDifference) {
1051               CircleFilterRowArray[RingRadiusInput][PointCount] = YInc;
1052               CircleFilterColumnArray[RingRadiusInput][PointCount] = XInc;
1053               PointCount++;
1054            }
1055         }
1056      }
1057
1058      /* number of points on given radius size */
1059      CircleFilterRowArray[RingRadiusInput][0] = PointCount-1;
1060      CircleFilterColumnArray[RingRadiusInput][0] = PointCount-1;
1061
1062   }
1063
1064   private static void MoatMaskCalc(double MoatMaskTempThreshold, double RingFitMaxRadiusDegree, int MoatSignCheckFlag) {
1065
1066      int ArrayInc,XInc,YInc;
1067      int NumberColumns=IRData_Remap_NumberColumns;
1068      int NumberRows=IRData_Remap_NumberRows;
1069    
1070      /* initialize arrays */
1071      for(YInc=0;YInc<NumberRows;YInc++) {
1072         for(XInc=0;XInc<NumberColumns;XInc++) {
1073            MoatMaskFlagField[YInc][XInc]=0;
1074            BlackWhiteFieldArray[YInc][XInc]=0;
1075         }
1076      }
1077  
1078      MeshGrid(MoatMaskTempThreshold,MoatSignCheckFlag);
1079      int NumberOfObjects = BWImage(8,NumberRows,NumberColumns);
1080    
1081      for(ArrayInc=1;ArrayInc<=NumberOfObjects;ArrayInc++) {
1082         double LatitudeMaximum = -90.0;
1083         double LatitudeMinimum = 90.0;
1084         double LongitudeMaximum = -180.0;
1085         double LongitudeMinimum = 180.0;
1086         for(YInc=0;YInc<NumberRows;YInc++) {
1087            for(XInc=0;XInc<NumberColumns;XInc++) {
1088               if(MoatMaskFlagField[YInc][XInc]==ArrayInc) {
1089                  double LatitudeValue = IRData_Remap_Latitude[YInc][XInc];
1090                  double LongitudeValue = IRData_Remap_Longitude[YInc][XInc];
1091                  if(LatitudeValue>LatitudeMaximum) {
1092                     LatitudeMaximum=LatitudeValue;
1093                  }
1094                  if(LongitudeValue>LongitudeMaximum) {
1095                     LongitudeMaximum=LongitudeValue;
1096                  }
1097                  if(LatitudeValue<LatitudeMinimum) {
1098                     LatitudeMinimum=LatitudeValue;
1099                  }
1100                  if(LongitudeValue<LongitudeMinimum) {
1101                     LongitudeMinimum=LongitudeValue;
1102                  }
1103               }
1104            }
1105         }
1106         double AverageLatitude = .5*(LatitudeMinimum+LatitudeMaximum);
1107         double FeatureLength = Math.sqrt(Math.pow((LatitudeMaximum-LatitudeMinimum),2)+
1108                                Math.pow((LongitudeMaximum-LongitudeMinimum)/Math.cos((PI*AverageLatitude)/180.0),2));
1109         int MoatFlagID = (FeatureLength>(RingFitMaxRadiusDegree*2.0)) ? 1 : 0;
1110         for(YInc=0;YInc<NumberRows;YInc++) {
1111            for(XInc=0;XInc<NumberColumns;XInc++) {
1112               if(MoatMaskFlagField[YInc][XInc]==ArrayInc) {
1113                  MoatMaskFlagField[YInc][XInc] = MoatFlagID;
1114               }
1115            }
1116         }
1117      }
1118
1119   }
1120
1121   private static void MeshGrid(double TemperatureThresholdValue, int MoatSignCheckFlag) {
1122
1123      int XInc,YInc;
1124
1125      for(XInc=0;XInc<IRData_Remap_NumberColumns;XInc++) {
1126         for(YInc=0;YInc<IRData_Remap_NumberRows;YInc++) {
1127            if(MoatSignCheckFlag==1) {
1128               BlackWhiteFieldArray[YInc][XInc] = (IRData_Remap_Temperature[YInc][XInc]> TemperatureThresholdValue) ? 1 : 0;
1129            } else {
1130               BlackWhiteFieldArray[YInc][XInc] = (IRData_Remap_Temperature[YInc][XInc]< TemperatureThresholdValue) ? 1 : 0;
1131            }
1132         }
1133      }
1134
1135   }
1136
1137   private static int BWImage(int ConnectednessValue, int NumberRows, int NumberColumns) {
1138      /* blackwhitefieldarray = BooleanValueArray; moatmaskflagfield = LabelImageArray */
1139    
1140      int IncVal,XInc,YInc;
1141      int B_Value,C_Value,D_Value,E_Value;
1142      int TableElements = 0;
1143      int LabelValue = 0;
1144      int[] LabelSetArray = new int[NumberRows*NumberColumns];
1145
1146      LabelSetArray[0] = 0;
1147    
1148      for(YInc=0;YInc<NumberRows;YInc++) {
1149         for(XInc=0;XInc<NumberColumns;XInc++) {
1150            if(BlackWhiteFieldArray[YInc][XInc]==1) {
1151               /* if A is an object */
1152               /* get the neighboring pixels B_Value, C_Value, D_Value, and E_Value */
1153               B_Value=0;
1154               C_Value=0;
1155               D_Value=0;
1156               E_Value=0;
1157               if(XInc>0) {
1158                  B_Value = Find(LabelSetArray,MoatMaskFlagField[YInc][XInc-1]);
1159               }
1160               if(YInc>0) {
1161                  C_Value = Find(LabelSetArray,MoatMaskFlagField[YInc-1][XInc]);
1162               }
1163               if((YInc>0)&&(XInc>0)) {
1164                  D_Value = Find(LabelSetArray,MoatMaskFlagField[YInc-1][XInc-1]);
1165               }
1166               if((YInc>0)&&(XInc>(NumberColumns-1))) {
1167                  E_Value = Find(LabelSetArray,MoatMaskFlagField[YInc-1][XInc+1]);
1168               }
1169               if(ConnectednessValue==4) {
1170                  /* apply 4 connectedness */
1171                  if((B_Value!=0)&&(C_Value!=0)) {
1172                     /* B_Value and C_Value are labeled */
1173                     if(B_Value==C_Value) {
1174                        MoatMaskFlagField[YInc][XInc] = B_Value;
1175                     } else {
1176                        LabelSetArray[C_Value] = B_Value;
1177                        MoatMaskFlagField[YInc][XInc] = B_Value;
1178                     }
1179                  } else if(B_Value!=0) {
1180                     /* B_Value is object but C_Value is not */
1181                     MoatMaskFlagField[YInc][XInc] = B_Value;
1182                  } else if(C_Value!=0) {
1183                     /* C_Value is object but B_Value is not */
1184                     MoatMaskFlagField[YInc][XInc] = C_Value;
1185                  } else {
1186                     /* B_Value, C_Value, D_Value not object - new object */
1187                     /* label and put into table */
1188                     TableElements++;
1189                     LabelSetArray[TableElements] = TableElements;
1190                     MoatMaskFlagField[YInc][XInc] = LabelSetArray[TableElements];
1191                  }
1192               } else if(ConnectednessValue==6) {
1193                  /* apply 6 connected ness */
1194                  if(D_Value!=0) {
1195                     /* D_Value object, copy label and move on */
1196                     MoatMaskFlagField[YInc][XInc] = D_Value;
1197                  } else if((B_Value!=0)&&(C_Value!=0)) {
1198                     /* B_Value and C_Value are labeled */
1199                     if(B_Value==C_Value) {
1200                        MoatMaskFlagField[YInc][XInc] = B_Value;
1201                     } else {
1202                        LabelValue = Math.min(B_Value,C_Value);
1203                        LabelSetArray[B_Value] = LabelValue;
1204                        LabelSetArray[C_Value] = LabelValue;
1205                        MoatMaskFlagField[YInc][XInc] = LabelValue;
1206                     }
1207                  } else if(B_Value!=0) {
1208                     /* B_Value is object but C_Value is not */
1209                     MoatMaskFlagField[YInc][XInc] = B_Value;
1210                  } else if(C_Value!=0) {
1211                     /* C_Value is object but B_Value is not */
1212                     MoatMaskFlagField[YInc][XInc] = C_Value;
1213                  } else {
1214                     /* B_Value, C_Value, D_Value not object - new object */
1215                     /* label and put into table */
1216                     TableElements++;
1217                     LabelSetArray[TableElements] = TableElements;
1218                     MoatMaskFlagField[YInc][XInc] = LabelSetArray[TableElements];
1219                  }
1220               } else if(ConnectednessValue==8) {
1221                  /* apply 8 connectedness */
1222                  if((B_Value!=0)||(C_Value!=0)||(D_Value!=0)||(E_Value!=0)) {
1223                     LabelValue = B_Value;
1224                     if(B_Value!=0) {
1225                        LabelValue = B_Value;
1226                     } else if(C_Value!=0) {
1227                        LabelValue = C_Value;
1228                     } else if(D_Value!=0) {
1229                        LabelValue = D_Value;
1230                     } else {
1231                        LabelValue = E_Value;
1232                     }
1233                     MoatMaskFlagField[YInc][XInc] = LabelValue;
1234                     if((B_Value!=0)&&(B_Value!=LabelValue)) {
1235                        LabelSetArray[B_Value] = LabelValue;
1236                     }
1237                     if((C_Value!=0)&&(C_Value!=LabelValue)) {
1238                        LabelSetArray[C_Value] = LabelValue;
1239                     }
1240                     if((D_Value!=0)&&(D_Value!=LabelValue)) {
1241                        LabelSetArray[D_Value] = LabelValue;
1242                     }
1243                     if((E_Value!=0)&&(E_Value!=LabelValue)) {
1244                        LabelSetArray[E_Value] = LabelValue;
1245                     }
1246                  } else {
1247                     /* label and put into table */
1248                     TableElements++;
1249                     LabelSetArray[TableElements] = TableElements;
1250                     MoatMaskFlagField[YInc][XInc] = LabelSetArray[TableElements];
1251                  }
1252               }
1253            } else {
1254               /* A is not an object so leave it */
1255               MoatMaskFlagField[YInc][XInc] = 0;
1256            }
1257         }
1258      }
1259    
1260      /* consolidate component table */
1261      for(IncVal=0;IncVal<=TableElements;IncVal++) {
1262         LabelSetArray[IncVal] = Find(LabelSetArray,IncVal);
1263      }
1264    
1265      /* run image through the look-up table */
1266      for(YInc=0;YInc<NumberRows;YInc++) {
1267         for(XInc=0;XInc<NumberColumns;XInc++) {
1268            MoatMaskFlagField[YInc][XInc] = LabelSetArray[MoatMaskFlagField[YInc][XInc]];
1269         }
1270      }
1271    
1272      /* count up the objects in the image */
1273      for(IncVal=0;IncVal<=TableElements;IncVal++) {
1274         LabelSetArray[IncVal] = 0;
1275      }
1276    
1277      for(YInc=0;YInc<NumberRows;YInc++) {
1278         for(XInc=0;XInc<NumberColumns;XInc++) {
1279            LabelSetArray[MoatMaskFlagField[YInc][XInc]]++;
1280         }
1281      }
1282    
1283      /* number the objects from 1 through n objects */
1284      int NumberObjects = 0;
1285      LabelSetArray[0] = 0;
1286      for(IncVal=1;IncVal<=TableElements;IncVal++) {
1287         if(LabelSetArray[IncVal]>0) {
1288            LabelSetArray[IncVal] = ++NumberObjects;
1289         }
1290      }
1291
1292      /* run through the look-up table again */
1293      for(YInc=0;YInc<NumberRows;YInc++) {
1294         for(XInc=0;XInc<NumberColumns;XInc++) {
1295            MoatMaskFlagField[YInc][XInc] = LabelSetArray[MoatMaskFlagField[YInc][XInc]];
1296         }
1297      }
1298
1299      return NumberObjects;
1300
1301   }
1302
1303   private static int Find(int InputArray[],int InputValue) {
1304      int IncVal = InputValue;
1305
1306      while (InputArray[IncVal]!=IncVal) {
1307         IncVal = InputArray[IncVal];
1308      }
1309
1310      return IncVal;
1311   }
1312   
1313   /**
1314    * This routine will determine the confidence scores for the spiral fitting
1315    * and ring fitting routines and calculate the best possible position for
1316    * the storm center based upon a combination of the two methods, if
1317    * available.
1318    *
1319    * If the ring fitting routine does not determine a good candidate position,
1320    * the spiral fitting routine will be used alone.  If the spiral fitting
1321    * routine candidate point is not "good", the forecast point will be
1322    * selected.
1323    *
1324    * @param FirstGuessLatitude First Guess latitude.
1325    * @param FirstGuessLongitude First Guess longitude.
1326    * @param SpiralCenterLatitude Spiral Analysis latitude at max location.
1327    * @param SpiralCenterLongitude Spiral Analysis longitude at max location.
1328    * @param SpiralCenterScoreValue Spiral Analysis Score value.
1329    * @param RingFitLatitude Ring Analysis latitude at max score location.
1330    * @param RingFitLongitude Ring Analysis longitude at max score location.
1331    * @param RingFitScoreValue Ring Analysis Score value.
1332    *
1333    * @return Array of four double values. The values represent: latitude of
1334    * final selected location, longitude of final selected location,
1335    * confidence score of final selected location, and the {@literal "method"}
1336    * used to determine the final selected location. Possible values for the
1337    * method: 1 for first guess, 2 for enhanced spiral analysis, and 5 for
1338    * combo ring/spiral analysis.
1339    */
1340   private static double[] CalcScores(double FirstGuessLatitude,
1341                                      double FirstGuessLongitude,
1342                                      double SpiralCenterLatitude,
1343                                      double SpiralCenterLongitude,
1344                                      double SpiralCenterScoreValue,
1345                                      double RingFitLatitude,
1346                                      double RingFitLongitude,
1347                                      double RingFitScoreValue)
1348   {
1349      int YInc;
1350      int FinalLocationSelectionMethodReturn = 0;
1351      int NumberOfSpirals = (int)SpiralCenterAnalysisField[0][0]+1;
1352
1353      double[] DistancePenaltyArray = new double[NumberOfSpirals];
1354      double[] InitialSpiralScoreArray = new double[NumberOfSpirals];
1355      double[] DistanceFromSpiralCenterArray = new double[NumberOfSpirals];
1356      double[] DistanceBonusArray = new double[NumberOfSpirals];
1357      double[] EnhancedSpiralScoreArray = new double[NumberOfSpirals];
1358      double[] FinalSpiralScoreArray = new double[NumberOfSpirals];
1359
1360      double FinalLatitudeReturn = -99.99;
1361      double FinalLongitudeReturn = -999.99;
1362      double FinalScoreReturn = 0.0;
1363      double DistanceValue;
1364      double MaximumForecastErrorDegree=1.0; /* max dist between fcst and final */
1365      double ExpectedMaximumForecastErrorDegree=MaximumForecastErrorDegree; /* expected error */
1366      
1367      /* max displacement */
1368      double MaximumAllowableDisplacement=ExpectedMaximumForecastErrorDegree*1.15;
1369      
1370      /* Spiral Center analysis weight */
1371      double SPIRALWEIGHT=10.0;
1372      
1373      /* RF distance penalty */
1374      double DISTPENALTYWEIGHT=(0.5/ExpectedMaximumForecastErrorDegree);
1375      
1376      /* Ring Fit bonus value */
1377      double PROXIMITYBONUS=4.5;
1378      
1379      /* RF bonus value threshold dist deg */
1380      double PROXIMITYTHRESH=0.25;
1381      
1382      /* combination score threshold value */
1383      double COMBOSCORETHRESH=15.0;
1384      
1385      /* convert distance in km to degrees */
1386      double KMperDegree=111.0;
1387      
1388      double SpiralCenterIndexMaximum = -999.99;
1389      double SpiralCenterMaximumLatitude = -99.99;
1390      double SpiralCenterMaximumLongitude = -999.99;
1391      double SpiralMaximumDistanceFromGuess = 999.99;
1392      double IntermediateRingScore = 0.0;
1393      double FinalSpiralScoreValue = 0.0;
1394      double MaximumSpiralScoreValue = -99.99;
1395      double MaximumSpiralScoreLatitude = -99.99;
1396      double MaximumSpiralScoreLongitude = -999.99;
1397
1398      /* Spiral Score Calculations */
1399      double LocalValue[] = Functions.distance_angle(FirstGuessLatitude,FirstGuessLongitude,
1400                            SpiralCenterLatitude,SpiralCenterLongitude,1);
1401      DistanceValue = LocalValue[0];
1402      double InitialSpiralScoreArrayScore = SpiralCenterScoreValue;
1403
1404      for(YInc=1;YInc<NumberOfSpirals;YInc++) {
1405         double LocalValue1[] = Functions.distance_angle(FirstGuessLatitude,FirstGuessLongitude,
1406                                                             SpiralCenterAnalysisField[YInc][1],
1407                                                             SpiralCenterAnalysisField[YInc][2],1);
1408         DistanceValue = LocalValue1[0];
1409         DistancePenaltyArray[YInc] = -DISTPENALTYWEIGHT*(DistanceValue/KMperDegree);
1410         InitialSpiralScoreArray[YInc] = SPIRALWEIGHT*(SpiralCenterAnalysisField[YInc][0]-InitialSpiralScoreArrayScore);
1411         double LocalValue2[] = Functions.distance_angle(SpiralCenterLatitude,SpiralCenterLongitude,
1412                                                             SpiralCenterAnalysisField[YInc][1],
1413                                                             SpiralCenterAnalysisField[YInc][2],1);
1414         DistanceValue = LocalValue2[0];
1415         DistanceFromSpiralCenterArray[YInc] = DistanceValue/KMperDegree;
1416         if(DistanceFromSpiralCenterArray[YInc]<=PROXIMITYTHRESH) {
1417            DistanceBonusArray[YInc] = PROXIMITYBONUS;
1418         } else {
1419            DistanceBonusArray[YInc] = 0.0;
1420         }
1421         EnhancedSpiralScoreArray[YInc] = InitialSpiralScoreArray[YInc]+DistancePenaltyArray[YInc];
1422         /*
1423          * System.out.printf("SpiralCenterAnalysisField :   YInc: %d  Lat: %f  Lon: %f Value: %f ",YInc,SpiralCenterAnalysisField[YInc][1],
1424          *                   SpiralCenterAnalysisField[YInc][2],SpiralCenterAnalysisField[YInc][0]);
1425          */
1426         /*
1427          * System.out.printf(" distPentaly=%f  DistanceFromSpiralCenterArray=%f DistanceBonusArray=%f   max=%f\n",
1428          *                   DistancePenaltyArray[YInc],DistanceFromSpiralCenterArray[YInc],
1429          *                   DistanceBonusArray[YInc],SpiralCenterIndexMaximum);
1430          */
1431         if(EnhancedSpiralScoreArray[YInc]>SpiralCenterIndexMaximum) {
1432            SpiralCenterIndexMaximum = EnhancedSpiralScoreArray[YInc];
1433            SpiralCenterMaximumLatitude = SpiralCenterAnalysisField[YInc][1];
1434            SpiralCenterMaximumLongitude = SpiralCenterAnalysisField[YInc][2];
1435         }
1436         FinalSpiralScoreArray[YInc] = InitialSpiralScoreArray[YInc]+DistancePenaltyArray[YInc]+DistanceBonusArray[YInc];
1437      }
1438      /*
1439       * System.out.printf("FirstGuessLatitude=%f  FirstGuessLongitude=%f SpiralCenterMaximumLatitude=%f SpiralCenterMaximumLongitude=%f\n",
1440       *                    FirstGuessLatitude,FirstGuessLongitude,SpiralCenterMaximumLatitude,SpiralCenterMaximumLongitude);
1441       */
1442      double LocalValue3[] = Functions.distance_angle(FirstGuessLatitude,FirstGuessLongitude,
1443                                                          SpiralCenterMaximumLatitude,SpiralCenterMaximumLongitude,1);
1444      DistanceValue = LocalValue3[0];
1445      SpiralMaximumDistanceFromGuess = DistanceValue/KMperDegree;
1446      /*
1447       * System.out.printf("SpiralMaximumDistanceFromGuess=%f MaximumAllowableDisplacement=%f\n",
1448       *                   SpiralMaximumDistanceFromGuess,MaximumAllowableDisplacement);
1449       */
1450      if(SpiralMaximumDistanceFromGuess<=MaximumAllowableDisplacement) {
1451         /* Ring Score Calculations */
1452         double LocalValue4[] = Functions.distance_angle(FirstGuessLatitude,FirstGuessLongitude,
1453                                                             RingFitLatitude,RingFitLongitude,1);
1454         DistanceValue = LocalValue4[0];
1455         for(YInc=1;YInc<NumberOfSpirals;YInc++) {
1456            double[] RingScoreReturn = FindRingScore(SpiralCenterAnalysisField[YInc][1],
1457                                                     SpiralCenterAnalysisField[YInc][2],
1458                                                     RingScoreAnalysisField);
1459            int RetErr = (int)RingScoreReturn[0];
1460            IntermediateRingScore = RingScoreReturn[1];
1461            if(RetErr==1) {
1462               FinalSpiralScoreValue = FinalSpiralScoreArray[YInc]+IntermediateRingScore;
1463               /*
1464                * System.out.printf("    FOUND  NumberOfSpirals=%d  lat=%f  lon=%f ringScore=%f  FinalSpiralScoreValue=%f MaximumSpiralScoreValue=%f\n",
1465                *                   YInc,SpiralCenterAnalysisField[YInc][1],SpiralCenterAnalysisField[YInc][2],
1466                *                   IntermediateRingScore,FinalSpiralScoreValue,MaximumSpiralScoreValue);
1467                */
1468               if(FinalSpiralScoreValue>MaximumSpiralScoreValue) {
1469                  MaximumSpiralScoreValue = FinalSpiralScoreValue;
1470                  MaximumSpiralScoreLatitude = SpiralCenterAnalysisField[YInc][1];
1471                  MaximumSpiralScoreLongitude = SpiralCenterAnalysisField[YInc][2];
1472               }
1473            }
1474         }
1475         /*
1476          * System.out.printf("MaximumSpiralScoreValue=%f RingPart=%f SpiralCenterIndexMaximum=%f\n",
1477          *                    MaximumSpiralScoreValue,IntermediateRingScore,SpiralCenterIndexMaximum);
1478          */
1479         if(MaximumSpiralScoreValue>=COMBOSCORETHRESH) {
1480            /* use Combo Method */
1481            FinalScoreReturn = MaximumSpiralScoreValue;
1482            FinalLatitudeReturn = MaximumSpiralScoreLatitude;
1483            FinalLongitudeReturn = MaximumSpiralScoreLongitude;
1484            FinalLocationSelectionMethodReturn = 5;
1485         }
1486         else {
1487            /* use ESP method */
1488            FinalScoreReturn = 1.0+SpiralCenterIndexMaximum;
1489            FinalLatitudeReturn = SpiralCenterMaximumLatitude;
1490            FinalLongitudeReturn = SpiralCenterMaximumLongitude;
1491            FinalLocationSelectionMethodReturn = 4;
1492         }
1493      }
1494      else {
1495         /* use Forecast Position */
1496         FinalScoreReturn = 0.0;
1497         FinalLatitudeReturn = FirstGuessLatitude;
1498         FinalLongitudeReturn = FirstGuessLongitude;
1499         FinalLocationSelectionMethodReturn = 1;
1500      }
1501
1502      return new double[] {
1503         FinalLatitudeReturn,
1504         FinalLongitudeReturn,
1505         FinalScoreReturn,
1506         (double)FinalLocationSelectionMethodReturn
1507      };
1508   }
1509   
1510   /**
1511    * Find Ring score at selected location (spiral analysis location)
1512    *
1513    * @param FirstGuessLatitude Latitude of search location.
1514    * @param FirstGuessLongitude Longitude of search location.
1515    * @param RingScoreAnalysisField - Array/Grid of Ring Analysis scores.
1516    *
1517    * @return Array of two doubles. The first value will be either -1 (if
1518    * nothing was found) or 1 (if found). The second value is the value at
1519    * search location in ring analysis grid (if the first value is 1).
1520    */
1521   private static double[] FindRingScore(double FirstGuessLatitude,
1522                                         double FirstGuessLongitude,
1523                                         double RingScoreAnalysisField[][])
1524   {
1525      int YInc = 1;
1526      int Ret_ID = -1;
1527      double MaxRingDistance = RingScoreAnalysisField[0][0];
1528      double RingFixScoreReturn = -99.9;
1529
1530      double LatitudeIncrement = Math.abs(IRData_Remap_Latitude[0][0]-IRData_Remap_Latitude[1][0])/2.0;
1531      double LongitudeIncrement = Math.abs(IRData_Remap_Longitude[0][1]-IRData_Remap_Longitude[0][0])/2.0;
1532      while(YInc<=MaxRingDistance) {
1533         if((Math.abs(RingScoreAnalysisField[YInc][1]-FirstGuessLatitude)<=LatitudeIncrement)&&
1534            (Math.abs(RingScoreAnalysisField[YInc][2]-FirstGuessLongitude)<=LongitudeIncrement)) {
1535            if(RingScoreAnalysisField[YInc][0]>RingFixScoreReturn) {
1536               RingFixScoreReturn = RingScoreAnalysisField[YInc][0];
1537            }
1538            Ret_ID = 1;
1539         }
1540         YInc++;
1541      }
1542
1543      return new double[] { (double)Ret_ID, RingFixScoreReturn };
1544   }
1545   
1546   /**
1547    * Determine method location scheme to use by examining various
1548    * empirically-defined confidence factors.
1549    *
1550    * Confidence factors will be derived,
1551    * with the "most confident" value used as the final automatically
1552    * determined storm position.
1553    *
1554    * @param ForecastLatitude NHC/JTWC interpolated latitude position.
1555    *
1556    * @param ForecastLongitude NHC/JTWC interpolated longitude position.
1557    *
1558    * @param RingSpiralLatitude Ring/Spiral Analysis latitude position.
1559    *
1560    * @param RingSpiralLongitude Ring/Spiral Analysis longitude position.
1561    *
1562    * @param RingSpiralScore Ring/Spiral Analysis confidence factor score.
1563    *
1564    * @param RingSpiralSelectionIDValue Ring/Spiral Analysis position
1565    *                                   derivation method.
1566    *
1567    * @return Array of three doubles. First value is latitude to be used,
1568    * second value is longitude to be used, and the third value is the
1569    * method used to determine storm position values.
1570    */
1571   private static double[] PickFinalLocation(int InputPositioningID,
1572                                             double ForecastLatitude,
1573                                             double ForecastLongitude,
1574                                             double RingSpiralLatitude,
1575                                             double RingSpiralLongitude,
1576                                             double RingSpiralScore,
1577                                             int RingSpiralSelectionIDValue)
1578   {
1579
1580      int PositioningMethodID = 0;
1581      int HistoryRecEyeScene;
1582      int HistoryRecCloudScene;
1583      double FinalLatitude = -99.99;
1584      double FinalLongitude = -99.99;
1585      double FinalScoreValue = 0.0;
1586      double HistoryRecFinalTno = 0.0;
1587      double MaximumHistoryRecCI = 0.0;
1588      double HistoryRecRawTno = 0.0;
1589      double HistoryRecCI = 0.0;
1590      boolean ForceAutoFixUseTF_forTesting = false;
1591      boolean FoundEyeSceneTF = false;
1592
1593      int CurDate = History.IRCurrentRecord.date;
1594      int CurTime = History.IRCurrentRecord.time;
1595      int CurCloudScene = History.IRCurrentRecord.cloudscene;
1596      int CurEyeScene = History.IRCurrentRecord.eyescene;
1597      double CurrentTime = Functions.calctime(CurDate,CurTime);
1598
1599      int NumRecsHistory = History.HistoryNumberOfRecords();
1600      double InitStrengthValue = Env.InitRawTValue;
1601      boolean LandFlagTF = Env.LandFlagTF;
1602
1603      int EyeSceneCount = 0;
1604
1605      if((Main.HistoryFileName!=null)||(ForceAutoFixUseTF_forTesting)) {
1606         HistoryRecFinalTno = 9.0;
1607         MaximumHistoryRecCI = 9.0;
1608         HistoryRecRawTno = 9.0;
1609      } else if(NumRecsHistory==0) {
1610         HistoryRecFinalTno = InitStrengthValue;
1611         MaximumHistoryRecCI = InitStrengthValue;
1612         HistoryRecRawTno = InitStrengthValue;
1613      } else {
1614         EyeSceneCount = 0;
1615         HistoryRecFinalTno = History.IRCurrentRecord.Traw;
1616         HistoryRecRawTno = History.IRCurrentRecord.Traw;
1617         MaximumHistoryRecCI = History.IRCurrentRecord.Traw;
1618         int XInc = 0;
1619         while(XInc<NumRecsHistory) {
1620            int RecDate = History.HistoryFile[XInc].date;
1621            int RecTime = History.HistoryFile[XInc].time;
1622            int RecLand = History.HistoryFile[XInc].land;
1623            double RecTnoRaw = History.HistoryFile[XInc].Traw;
1624            double HistoryRecTime = Functions.calctime(RecDate,RecTime);
1625            boolean LandCheckTF = true;
1626            if(((LandFlagTF)&&(RecLand==1))||(RecTnoRaw<1.0)) {
1627               LandCheckTF = false;
1628            }
1629            if((HistoryRecTime<CurrentTime)&&(LandCheckTF)) {
1630               HistoryRecCI = History.HistoryFile[XInc].CI;
1631               HistoryRecFinalTno = History.HistoryFile[XInc].Tfinal;
1632               HistoryRecRawTno = History.HistoryFile[XInc].Traw;
1633               HistoryRecEyeScene = History.HistoryFile[XInc].eyescene;
1634               HistoryRecCloudScene = History.HistoryFile[XInc].cloudscene;
1635               if(HistoryRecCI>MaximumHistoryRecCI) {
1636                  MaximumHistoryRecCI = HistoryRecCI;
1637               }
1638               if((HistoryRecEyeScene<3)||((HistoryRecCloudScene==1)&&(HistoryRecEyeScene==3))) {
1639                  EyeSceneCount = EyeSceneCount+1;
1640                  /* checking for eye or embedded center scene types */
1641                  if(EyeSceneCount==3) {
1642                     FoundEyeSceneTF = true;
1643                  }
1644               } else {
1645                  EyeSceneCount = 0;
1646               }
1647            }
1648            XInc++;
1649         }
1650         if(((CurEyeScene<3)||((CurCloudScene==1)&&(CurEyeScene==3)))&&(EyeSceneCount==2)) {
1651            FoundEyeSceneTF = true;
1652         }
1653      }
1654
1655      /* System.out.printf("%d %d : ",CurDate,CurTime); */
1656      /* System.out.printf("RingSpiralScore=%f HistoryRecFinalTno=%f\n",RingSpiralScore,HistoryRecFinalTno); */
1657      /* System.out.printf("HistoryRecRawTno=%f MaximumHistoryRecCI=%f\n",HistoryRecRawTno,MaximumHistoryRecCI); */
1658    
1659      /*
1660       * check score for developing systems
1661       * (MaximumHistoryRecCI<5.0) starting at 3.0 or
1662       * check score for weakeining systems
1663       * (MaximumHistoryRecCI>=5.0) only above 3.5
1664       */
1665      if(((HistoryRecRawTno>=3.0)&&(MaximumHistoryRecCI<5.0))||
1666         ((HistoryRecRawTno>=3.5)&&(MaximumHistoryRecCI>=5.0))) {
1667    
1668         if((HistoryRecFinalTno<=4.5)&&(RingSpiralScore<1.0)&&
1669            (!FoundEyeSceneTF)) RingSpiralScore = -99.99;
1670
1671         /* System.out.printf(" FINALRingSpiralScore=%f ",RingSpiralScore); */
1672    
1673         if(RingSpiralScore>0.0) {
1674            /* use Spiral/Ring methodology for center point */
1675            FinalLatitude = RingSpiralLatitude;
1676            FinalLongitude = RingSpiralLongitude;
1677            FinalScoreValue = RingSpiralScore;
1678            PositioningMethodID = RingSpiralSelectionIDValue;
1679         } else {
1680            /* CDO... can't find anything to focus on */
1681            FinalLatitude = ForecastLatitude;
1682            FinalLongitude = ForecastLongitude;
1683            FinalScoreValue = 0.0;
1684            PositioningMethodID = InputPositioningID;
1685         }
1686      } else {
1687         /*
1688          * current Tfinal is less than 3.5 or current scene
1689          * is not an eye or embedded center
1690          * WILL USE FORECAST POSITION FOR AODT ANALYSIS
1691          */
1692         FinalLatitude = ForecastLatitude;
1693         FinalLongitude = ForecastLongitude;
1694         FinalScoreValue = 0.0;
1695         PositioningMethodID = InputPositioningID;
1696      }
1697    
1698      return new double[] { FinalLatitude, FinalLongitude, FinalScoreValue, (double)PositioningMethodID};
1699
1700   }
1701 
1702
1703   
1704   /**
1705    * Calls routines to setup transformation, transform, data move.
1706    *
1707    * Input data provided with global variable arrays containing original
1708    * and transformed arrays.
1709    */
1710   private static void RemapData() {
1711      /*
1712       * The following routines were originally developed by Dave Santek
1713       * of UW/SSEC and were added to the ADT under permission.
1714       * If executed, an array of latitude and longitude position arrays
1715       * will be remapped to a rectilinear projection for Tony Wimmers routines
1716       */
1717
1718      int NumberOfCorners = 0;
1719      int LineSplineValue = 3;
1720      int ElementSplineValue = LineSplineValue;
1721    
1722      tiff_vars.in_elems = Data.IRData_NumberColumns;
1723      tiff_vars.in_lines = Data.IRData_NumberRows;
1724    
1725      /* System.out.printf("elems=%d lines=%d\n", tiff_vars.in_elems,tiff_vars.in_lines); */
1726
1727      /* Uinit( ); */
1728      /*
1729      dis_vars.xrectl = (double)tiff_vars.in_lines;
1730      dis_vars.xrecte = (double)tiff_vars.in_elems;
1731      System.out.printf("xrectl=%f xrecte=%f\n",dis_vars.xrectl,dis_vars.xrecte);
1732      */
1733      
1734      DetermineDest( );
1735      /* System.out.printf("after determineDest\n"); */
1736      
1737      NumberOfCorners = Init(LineSplineValue,ElementSplineValue);
1738
1739      System.out.printf("number of corners=%d\n",NumberOfCorners);
1740      Corner(NumberOfCorners,LineSplineValue,ElementSplineValue);
1741      /* System.out.printf("after corners\n"); */
1742
1743      if((ElementSplineValue>1)||(LineSplineValue>1)) {
1744         DoMap(NumberOfCorners,LineSplineValue,ElementSplineValue);
1745      }
1746
1747   }
1748   
1749   /**
1750    * Interpolate between two arrays of different size.
1751    */
1752   private static void DetermineDest() {
1753   int XInc,YInc;
1754      int XSizeMax= Data.IRData_NumberColumns-1;
1755      int YSizeMax= Data.IRData_NumberRows-1;
1756       
1757      double NWCornerLatitude = Data.IRData_Latitude[0][0];
1758      double NWCornerLongitude = Data.IRData_Longitude[0][0];
1759      double NECornerLatitude = Data.IRData_Latitude[0][XSizeMax];
1760      double NECornerLongitude = Data.IRData_Longitude[0][XSizeMax];
1761      double SWCornerLatitude = Data.IRData_Latitude[YSizeMax][0];
1762      double SWCornerLongitude = Data.IRData_Longitude[YSizeMax][0];
1763      double SECornerLatitude = Data.IRData_Latitude[YSizeMax][XSizeMax];
1764      double SECornerLongitude = Data.IRData_Longitude[YSizeMax][XSizeMax];
1765
1766      /* crosses dateline check */
1767      if((NWCornerLongitude<NECornerLongitude)||(SWCornerLongitude<SECornerLongitude)) {
1768         NWCornerLongitude = NWCornerLongitude+360.0;
1769         if(SWCornerLongitude<SECornerLongitude) {
1770            SWCornerLongitude = SWCornerLongitude+360.0;
1771         }
1772         /* System.out.printf("DATELINE CROSS\n"); */
1773         for(XInc=0;XInc<XSizeMax;XInc++) {
1774            for(YInc=0;YInc<YSizeMax;YInc++) {
1775               double DataValue = Data.IRData_Longitude[YInc][XInc];
1776               if(DataValue<0.0) {
1777                 Data.IRData_Longitude[YInc][XInc] = (float)(DataValue+360.0);
1778               }
1779            }
1780         }
1781         NWCornerLatitude = Data.IRData_Latitude[0][0];
1782         NWCornerLongitude = Data.IRData_Longitude[0][0];
1783         NECornerLatitude = Data.IRData_Latitude[0][XSizeMax];
1784         NECornerLongitude = Data.IRData_Longitude[0][XSizeMax];
1785         SWCornerLatitude = Data.IRData_Latitude[YSizeMax][0];
1786         SWCornerLongitude = Data.IRData_Longitude[YSizeMax][0];
1787         SECornerLatitude = Data.IRData_Latitude[YSizeMax][XSizeMax];
1788         SECornerLongitude = Data.IRData_Longitude[YSizeMax][XSizeMax];
1789      }
1790
1791      double MaximumLatitudeValue = Math.min(NWCornerLatitude,NECornerLatitude);
1792      double MinimumLatitudeValue = Math.max(SWCornerLatitude,SECornerLatitude);
1793      double MaximumLongitudeValue = Math.max(NWCornerLongitude,SWCornerLongitude);
1794      double MinimumLongitudeValue = Math.min(NECornerLongitude,SECornerLongitude);
1795
1796      double LatitudeIncrement = (MaximumLatitudeValue-MinimumLatitudeValue)/
1797                                 (double) Data.IRData_NumberColumns;
1798      double LongitudeIncrement = (MaximumLongitudeValue-MinimumLongitudeValue)/
1799                                  (double) Data.IRData_NumberRows;
1800      double MaximumIncrementValue = Math.max(LatitudeIncrement,LongitudeIncrement);
1801        System.out.printf("REMAPPING INFO\n");
1802        System.out.printf("Source Array Bounds\n");
1803        System.out.printf("  NW Corner : %7.2f/%7.2f\n",NWCornerLatitude,NWCornerLongitude);
1804        System.out.printf("  NE Corner : %7.2f/%7.2f\n",NECornerLatitude,NECornerLongitude);
1805        System.out.printf("  SW Corner : %7.2f/%7.2f\n",SWCornerLatitude,SWCornerLongitude);
1806        System.out.printf("  SE Corner : %7.2f/%7.2f\n",SECornerLatitude,SECornerLongitude);
1807        System.out.printf("Destination Array Bounds\n");
1808        System.out.printf("  Max Lat/Lon: %7.2f/%7.2f\n",MaximumLatitudeValue,MaximumLongitudeValue);
1809        System.out.printf("  Min Lat/Lon: %7.2f/%7.2f\n",MinimumLatitudeValue,MinimumLongitudeValue);
1810        System.out.printf("  Inc Lat/Lon:   %5.3f/  %5.3f\n",MaximumIncrementValue,MaximumIncrementValue);
1811       
1812      tiff_vars.out_lines = (int)((MaximumLatitudeValue-MinimumLatitudeValue)/MaximumIncrementValue);
1813      tiff_vars.out_elems = (int)((MaximumLongitudeValue-MinimumLongitudeValue)/MaximumIncrementValue);
1814      for (YInc = 0; YInc < tiff_vars.out_lines; YInc++) {
1815         double LatitudeValue = MaximumLatitudeValue-(YInc*MaximumIncrementValue);
1816         for (XInc = 0; XInc < tiff_vars.out_elems; XInc++) {
1817            double LongitudeValue = MaximumLongitudeValue-(XInc*MaximumIncrementValue);
1818            IRData_Remap_Latitude[YInc][XInc] = LatitudeValue;
1819            IRData_Remap_Longitude[YInc][XInc] = LongitudeValue;
1820         }
1821      }
1822
1823      IRData_Remap_NumberColumns = tiff_vars.out_elems;
1824      IRData_Remap_NumberRows = tiff_vars.out_lines;
1825   }
1826   
1827   /**
1828    * Compute number of corners for transformation and block sizes.
1829    *
1830    * @param LineSplineInput Spline function for line values.
1831    * @param ElementSplineInput Spline function for element values.
1832    *
1833    * @return Total number of corners to interpolate.
1834    */
1835   private static int Init(int LineSplineInput, int ElementSplineInput) {
1836      int NumberOfCorners = 0;
1837    
1838      remap_vars.nspl = (tiff_vars.out_elems+ElementSplineInput-1)/ElementSplineInput;
1839      remap_vars.nspe = (tiff_vars.out_lines+LineSplineInput-1)/LineSplineInput;
1840    
1841      remap_vars.ncl = remap_vars.nspl+1;
1842      remap_vars.nce = remap_vars.nspe+1;
1843    
1844      if(((tiff_vars.out_elems+ElementSplineInput-1)%ElementSplineInput)==0) {
1845         remap_vars.ncl = remap_vars.nspl;
1846      }
1847      if(((tiff_vars.out_lines+LineSplineInput-1)%LineSplineInput)==0) {
1848         remap_vars.nce = remap_vars.nspe;
1849      }
1850    
1851      NumberOfCorners = remap_vars.ncl*remap_vars.nce;
1852    
1853      remap_vars.in_bfw = Math.max(MINBFW,MINBLKSIZ*tiff_vars.in_elems);
1854      remap_vars.out_bfw = Math.max(MINBFW,Math.max(LineSplineInput,MINBLKSIZ)*tiff_vars.out_elems);
1855    
1856      remap_vars.slb = remap_vars.in_bfw/tiff_vars.in_elems;
1857      remap_vars.dlb = ((remap_vars.out_bfw/tiff_vars.out_elems)/LineSplineInput)*LineSplineInput;
1858
1859      return NumberOfCorners;
1860   }
1861   
1862   /**
1863    * Compute transformations at corners.
1864    *
1865    * Operates on the {@link #LineCoordinateArray} and
1866    * {@link #ElementCoordinateArray}.
1867    *
1868    * @param NumberOfCornersInput Total number of corners to interpolate.
1869    * @param LineSplineInput Spline function for line values.
1870    * @param ElementSplineInput Spline function for element values.
1871    */
1872   private static void Corner(int NumberOfCornersInput,
1873                              int LineSplineInput,
1874                              int ElementSplineInput)
1875   {
1876      int XInc;
1877      int YInc;
1878      int ArrayInc;
1879
1880      /* initialize array of corners */
1881      for (ArrayInc=0;ArrayInc<NumberOfCornersInput;ArrayInc++) {
1882         LineCoordinateArray[ArrayInc] = -99.9;
1883         ElementCoordinateArray[ArrayInc] = -99.9;
1884      }
1885    
1886      /* loop through destination file and record source coords */
1887      ArrayInc = -1;
1888      int NumberLines = tiff_vars.out_lines+LineSplineInput-1;
1889      int NumberElements = tiff_vars.out_elems+ElementSplineInput-1;
1890    
1891      for (YInc=0;YInc<NumberLines;YInc=YInc+LineSplineInput) {
1892         int LineValue = YInc;
1893         for (XInc=0;XInc<NumberElements;XInc=XInc+ElementSplineInput) {
1894            /* System.out.printf("yinc=%d xinc=%d... \n",YInc,XInc); */
1895            int ElementValue = XInc;
1896            int[] UMapReturn = UMap(LineValue,ElementValue);
1897            int RetErr = UMapReturn[0];
1898            int SplineLineValue = UMapReturn[1];
1899            int SplineElementValue = UMapReturn[2];
1900            /* System.out.printf("spline line=%d element=%d \n ",SplineLineValue,SplineElementValue); */
1901            if((ElementSplineInput==1)&&(LineSplineInput==1)) {
1902               IRData_Remap_Temperature[LineValue][ElementValue]=
1903                       Data.IRData_Temperature[SplineLineValue][SplineElementValue];
1904            } else {
1905               ++ArrayInc;
1906               if(RetErr==0) {
1907                  LineCoordinateArray[ArrayInc] = SplineLineValue;
1908                  ElementCoordinateArray[ArrayInc] = SplineElementValue;
1909               }
1910            }
1911         }
1912      }
1913
1914   }
1915   
1916   /**
1917    * Provide coordinates between original point and transformed point.
1918    *
1919    * @param LineValueInput Original line coordinate.
1920    * @param ElementValueInput Ooriginal element coordinate.
1921    *
1922    * @return Array containing three values: possible error code, interpolated
1923    * line coordinate, and interpolated element coordinate.
1924    */
1925   private static int[] UMap(int LineValueInput, int ElementValueInput) {
1926
1927      /* dest array line position value */
1928      double DestinationLatitude;
1929
1930      /* dest array element position value */
1931      double DestinationLongitude;
1932    
1933      /*
1934       * Convert destination LineValueInput, ElementValueInput to
1935       * source coords LineValue_Return, ElementValue_Return
1936       */
1937      LineValueInput = Math.min(LineValueInput,tiff_vars.out_lines-1);
1938      ElementValueInput = Math.min(ElementValueInput,tiff_vars.out_elems-1);
1939      DestinationLatitude = IRData_Remap_Latitude[LineValueInput][ElementValueInput];
1940      DestinationLongitude = IRData_Remap_Longitude[LineValueInput][ElementValueInput];
1941    
1942      int[] FindPointReturn = FindPoint(DestinationLatitude,DestinationLongitude);
1943      int RetErr = FindPointReturn[0];
1944      int LineValueReturn = FindPointReturn[1];
1945      int ElementValueReturn = FindPointReturn[2];
1946
1947      return new int[] { RetErr, LineValueReturn, ElementValueReturn };
1948   }
1949   
1950   /**
1951    * Find specific lat/lon location in array and return index values.
1952    *
1953    * @param latitude Latitude value.
1954    * @param longitude Longitude value.
1955    *
1956    * @return Array of three values. The first value is the status (-1 for
1957    * error, 0 for ok), the second value is array line value of lat/lon
1958    * input, and the third value is array element value of lat/lon input.
1959    */
1960   private static int[] FindPoint(double latitude, double longitude) {
1961      int RetErr;
1962      int XInc;
1963      int IndexValue;
1964      int YLineValueReturn;
1965      int XElementValueReturn;
1966
1967      double[] CornerLatitudeArray = new double[4];
1968      double[] CornerLongitudeArray = new double[4];
1969      double[] NSEWDistanceValuesArray = new double[4];
1970
1971      /* found point logical */
1972      boolean FoundPointTF = false;
1973      
1974      /* out of bounds flag logical */
1975      boolean OutOfBoundsTF = false;
1976      
1977      /* found latitude value logical */
1978      boolean FoundLatitudeTF = false;
1979      
1980      /* found longitude value logical */
1981      boolean FoundLongitudeTF = false;
1982      
1983
1984      int NumberElements = tiff_vars.in_elems;
1985      int NumberLines = tiff_vars.in_lines;
1986      int ElementValue = 0;
1987      int LineValue = 0;
1988      double PreviousDistance = 9999.9;
1989
1990      for(XInc=0;XInc<4;XInc++) {
1991         NSEWDistanceValuesArray[XInc] = 0.0;
1992         CornerLatitudeArray[XInc] = 0.0;
1993         CornerLongitudeArray[XInc] = 0.0;
1994      }
1995      while((!FoundPointTF)&&(!OutOfBoundsTF)) {
1996         CornerLatitudeArray[0] = Data.IRData_Latitude[LineValue][ElementValue];
1997         CornerLongitudeArray[0] = Data.IRData_Longitude[LineValue][ElementValue];
1998         CornerLatitudeArray[3] = Data.IRData_Latitude[LineValue+1][ElementValue+1];
1999         CornerLongitudeArray[3] = Data.IRData_Longitude[LineValue+1][ElementValue+1];
2000         /*
2001          * System.out.printf("x=%d  y=%d  : CornerLatitudeArray0=%f CornerLongitudeArray0=%f ",
2002          *                   " CornerLatitudeArray3=%f CornerLongitudeArray3=%f\n",
2003          *                   ElementValue,LineValue,
2004          *                   CornerLatitudeArray[0],CornerLongitudeArray[0],
2005          *                   CornerLatitudeArray[3],CornerLongitudeArray[3]);
2006          */
2007         if((longitude>CornerLongitudeArray[0])||
2008            (longitude<CornerLongitudeArray[3])) {
2009            FoundLongitudeTF = false;
2010            if(longitude<CornerLongitudeArray[3]) {
2011               ElementValue++;
2012            } else {
2013               if(longitude>CornerLongitudeArray[0]) {
2014                  ElementValue--;
2015               }
2016            }
2017         } else {
2018            FoundLongitudeTF = true;
2019         }
2020         if((latitude>CornerLatitudeArray[0])||
2021            (latitude<CornerLatitudeArray[3])) {
2022            FoundLatitudeTF = false;
2023            if(latitude<CornerLatitudeArray[3]) {
2024               LineValue++;
2025            } else {
2026               if(latitude>CornerLatitudeArray[0]) {
2027                  LineValue--;
2028               }
2029            }
2030         } else {
2031            FoundLatitudeTF = true;
2032         }
2033         double LocalValue1[] = Functions.distance_angle(latitude,longitude,
2034                                    CornerLatitudeArray[0],CornerLongitudeArray[0],1);
2035         double DistanceValue = LocalValue1[0];
2036
2037         /*
2038         ** System.out.printf("distance : latitude=%f longitude=%f",
2039         **                   " CornerLatitudeArray0=%f CornerLongitudeArray0=%f",
2040         **                   " dist=%f angle=%f\n",latitude,longitude,
2041         **                    CornerLatitudeArray[0],CornerLongitudeArray[0],
2042         **                    DistanceValue,AngleValue);
2043         */
2044         if(FoundLatitudeTF&&FoundLongitudeTF) {
2045            FoundPointTF = true;
2046         }
2047         if(PreviousDistance<=DistanceValue) {
2048            FoundPointTF = true;
2049         }
2050         if((ElementValue<0)||
2051             (ElementValue>(NumberElements-2))) {
2052            OutOfBoundsTF = true;
2053         }
2054         if((LineValue<0)||(LineValue>(NumberLines-2))) {
2055            OutOfBoundsTF = true;
2056         }
2057         PreviousDistance = DistanceValue;
2058      }
2059      if(FoundPointTF) {
2060         CornerLatitudeArray[1] = Data.IRData_Latitude[LineValue][ElementValue+1];
2061         CornerLongitudeArray[1] = Data.IRData_Longitude[LineValue][ElementValue+1];
2062         CornerLatitudeArray[2] = Data.IRData_Latitude[LineValue+1][ElementValue];
2063         CornerLongitudeArray[2] = Data.IRData_Longitude[LineValue+1][ElementValue];
2064
2065         double LocalValue2[] = Functions.distance_angle(latitude,longitude,
2066                                    CornerLatitudeArray[0],CornerLongitudeArray[0],1);
2067         NSEWDistanceValuesArray[0] = LocalValue2[0];
2068         double LocalValue3[] = Functions.distance_angle(latitude,longitude,
2069                                    CornerLatitudeArray[1],CornerLongitudeArray[1],1);
2070         NSEWDistanceValuesArray[1] = LocalValue3[0];
2071         IndexValue = (NSEWDistanceValuesArray[0]<NSEWDistanceValuesArray[1]) ? 0 : 1;
2072         double LocalValue4[] = Functions.distance_angle(latitude,longitude,
2073                                    CornerLatitudeArray[2],CornerLongitudeArray[2],1);
2074         NSEWDistanceValuesArray[2] = LocalValue4[0];
2075         IndexValue = (NSEWDistanceValuesArray[IndexValue]<NSEWDistanceValuesArray[2]) ? IndexValue : 2;
2076         double LocalValue5[] = Functions.distance_angle(latitude,longitude,
2077                                    CornerLatitudeArray[3],CornerLongitudeArray[3],1);
2078         NSEWDistanceValuesArray[3] = LocalValue5[0];
2079         IndexValue = (NSEWDistanceValuesArray[IndexValue]<NSEWDistanceValuesArray[3]) ? IndexValue : 3;
2080
2081         XElementValueReturn = ((IndexValue==0)||(IndexValue==2)) ? ElementValue : ElementValue+1;
2082         YLineValueReturn = ((IndexValue==0)||(IndexValue==1)) ? LineValue : LineValue+1;
2083         RetErr = 0;
2084      } else {
2085         XElementValueReturn = -1;
2086         YLineValueReturn = -1;
2087         RetErr = -1;
2088      }
2089
2090      return new int[] { RetErr, YLineValueReturn, XElementValueReturn };
2091   }
2092
2093   private static int DoMap(int NumberOfCornersInput, int LineSplineInput, int ElementSplineInput) {
2094      /* LineCoordsArrayInput = LineCoordinateArray */
2095      /* ElementCoordsArrayInput = ElementCoordinateArray */
2096
2097      /* line increment */
2098      int[] LineIncArray1 = { -2, -2, 0, 2, 2, 2, 0, -2 };
2099      
2100      /* line increment */
2101      int[] LineIncArray2 = { -1, -1, 0, 1, 1, 1, 0, -1 };
2102      
2103      /* ele increment*/
2104      int[] ElementIncArray1 = { 0, 2, 2, 2, 0, -2, -2, -2 };
2105      
2106      /* ele increment*/
2107      int[] ElementIncArray2 = { 0, 1, 1, 1, 0, -1, -1, -1 };
2108      
2109   
2110      int BufferSize = tiff_vars.in_lines*tiff_vars.in_elems;
2111   
2112      double[] SourceArray = new double[BufferSize];
2113      double[] InterpArray = new double[BufferSize];
2114      double[] TempCharArray = new double[NumberOfCornersInput];
2115
2116      int Ret_ID;
2117      int XInc,YInc,IncVal;
2118      int ArrayIndexValue;
2119      int ArrayPointValue;
2120      int SplinePixelValue;
2121      int SourceBlockValue;
2122      int SourceBlockInc;
2123      int SplineInc;
2124      int OffsetValue;
2125      int OffsetValue0;
2126      int OffsetValue1;
2127      int OffsetValue2;
2128      int OffsetValueL = 0;
2129      int OffsetValueE = 0;
2130      int CornerPointer;
2131      int AccumulatedLines;
2132      int FirstCornerPointer;
2133      int FirstCornerPointerHold;
2134      int LastCornerPointer;
2135      int ValX;
2136      int LineMinimumFinal;
2137      int LineMaximumFinal;
2138      int ElementMinimumFinal;
2139      int ElementMaximumFinal;
2140      int EdgeFixID;
2141      int SourceLine = 0;
2142      int SourceElement = 0;
2143      int SourceOffsetValue;
2144      int MaximumSourceLine;
2145
2146      double LineValueA;
2147      double LineValueB;
2148      double LineValueC;
2149      double ElementValueA;
2150      double ElementValueB;
2151      double ElementValueC;
2152      double LineValueACounter;
2153      double LineValueBCounter;
2154      double LineValueCCounter;
2155      double ElementValueACounter;
2156      double ElementValueBCounter;
2157      double ElementValueCCounter;
2158      double LineValueACounter0;
2159      double ElementValueACounter0;
2160      double LineULValue;
2161      double LineURValue;
2162      double LineLLValue;
2163      double LineLRValue;
2164      double LineMinimum;
2165      double LineMaximum;
2166      double ElementULValue;
2167      double ElementURValue;
2168      double ElementLLValue;
2169      double ElementLRValue;
2170      double ElementMinimum;
2171      double ElementMaximum;
2172      boolean ReadAllTF = true;
2173      boolean SKIP;
2174      boolean LABEL30;
2175      boolean LABEL38;
2176
2177      if((SourceArray==null)||(InterpArray==null)||(TempCharArray==null)) {
2178         Ret_ID = -1;
2179      } else {
2180         SplinePixelValue = LineSplineInput*ElementSplineInput;
2181    
2182         for(YInc=0;YInc<tiff_vars.in_lines;YInc++) {
2183            for(XInc=0;XInc<tiff_vars.in_elems;XInc++) {
2184               ArrayIndexValue = (YInc*tiff_vars.in_elems)+XInc;
2185               SourceArray[ArrayIndexValue] = Data.IRData_Temperature[YInc][XInc];
2186            }
2187         }
2188    
2189         /* System.out.printf("remap nce=%d ncl=%d\n",remap_vars.nce,remap_vars.ncl); */
2190         for(YInc=1;YInc<remap_vars.nce+1;YInc++) {
2191            for(XInc=1;XInc<remap_vars.ncl+1;XInc++) {
2192               ArrayPointValue = IND(YInc,XInc);
2193
2194               if(LineCoordinateArray[ArrayPointValue-1]==(double)-99.9) {
2195                  double LineSum = 0.0;
2196                  double ElementSum = 0.0;
2197                  int PointCounter = 0;
2198                  for(IncVal=0;IncVal<8;IncVal++) {
2199                     SKIP = false;
2200                     int YInc1 = YInc+LineIncArray1[IncVal];
2201                     int YInc2 = YInc+LineIncArray2[IncVal];
2202                     if((YInc1<1)||(YInc2<1)) {
2203                        SKIP = true;
2204                     }
2205                     if((YInc1>remap_vars.nce)||(YInc2>remap_vars.nce)) {
2206                        SKIP = true;
2207                     }
2208    
2209                     int XInc1=XInc+ElementIncArray1[IncVal];
2210                     int XInc2=XInc+ElementIncArray2[IncVal];
2211
2212                     if((XInc1<1)||(XInc2<1)) {
2213                        SKIP = true;
2214                     }
2215                     if((XInc1>remap_vars.ncl)||(XInc2>remap_vars.ncl)) {
2216                        SKIP = true;
2217                     }
2218    
2219                     int ValueIND1 = IND(YInc1,XInc1);
2220                     int ValueIND2 = IND(YInc2,XInc2);
2221                     
2222                     if(!SKIP) {
2223                        if((LineCoordinateArray[ValueIND1-1]==-99.9)||(LineCoordinateArray[ValueIND2-1]==-99.9)) { 
2224                           SKIP = true;
2225                        }
2226                        if(TempCharArray[ValueIND1-1]!= 0) {
2227                           SKIP = true;
2228                        }
2229                        if(TempCharArray[ValueIND2-1]!= 0) {
2230                           SKIP = true;
2231                        }
2232                        if(!SKIP) {
2233                           PointCounter = PointCounter+1;
2234                           LineSum = LineSum+
2235                                  ((double)(2.0*LineCoordinateArray[ValueIND2-1]))-LineCoordinateArray[ValueIND1-1];
2236                           ElementSum = ElementSum+
2237                                  ((double)2.0*ElementCoordinateArray[ValueIND2-1])-ElementCoordinateArray[ValueIND1-1];
2238                        } /* SKIP:; */
2239                     } /* SKIP:; */
2240                  }
2241                      
2242                  if(PointCounter>0) {
2243                     LineCoordinateArray[ArrayPointValue-1] = LineSum/PointCounter;
2244                     ElementCoordinateArray[ArrayPointValue-1] = ElementSum/PointCounter;
2245                     TempCharArray[ArrayPointValue-1] = 1;
2246                  }
2247               }
2248            }
2249         }
2250             
2251         /* Loop through by destination blocks */
2252         int BlockCounter = 0;
2253    
2254         for(IncVal=1;IncVal<tiff_vars.out_lines+1;IncVal=IncVal+remap_vars.dlb) {
2255            /* Accumulated lines/block */
2256            AccumulatedLines = BlockCounter*remap_vars.dlb;
2257    
2258            /* Pointer to first corner of splines for this dest block */
2259            FirstCornerPointer = (BlockCounter*remap_vars.ncl*remap_vars.dlb)/LineSplineInput;
2260            FirstCornerPointerHold = FirstCornerPointer;
2261    
2262            /* Pointer to last corner for this dest block */
2263            LastCornerPointer = (((BlockCounter+1)*remap_vars.ncl*remap_vars.dlb)/LineSplineInput)-1;
2264            LastCornerPointer=Math.min(LastCornerPointer,NumberOfCornersInput-remap_vars.ncl);
2265    
2266            /* For each destination block loop through entire source */
2267            for(SourceBlockInc=1;SourceBlockInc<tiff_vars.in_lines+1;SourceBlockInc=SourceBlockInc+remap_vars.slb) {
2268               MaximumSourceLine = Math.min(tiff_vars.in_lines,(SourceBlockInc+remap_vars.slb-1));
2269               ReadAllTF = false;
2270    
2271               /* Loop through splines and move any data */
2272    
2273               FirstCornerPointer = FirstCornerPointerHold;
2274               SourceBlockValue = 0;
2275               while(FirstCornerPointer<LastCornerPointer) {
2276                  for(SplineInc=1;SplineInc<remap_vars.nspl+1;SplineInc++) {
2277                     LABEL30 = false;
2278                     OffsetValue0 = (SourceBlockValue/remap_vars.nspl)*LineSplineInput*tiff_vars.out_elems;
2279                     OffsetValue1 = (SourceBlockValue%remap_vars.nspl)*ElementSplineInput;
2280                     OffsetValue2 = OffsetValue1+OffsetValue0+1;
2281                     CornerPointer = FirstCornerPointer+SplineInc-1;
2282    
2283                     /* Get 4 corners in line space and check for out of bounds */
2284    
2285                     ValX = remap_vars.ncl;
2286                     LineULValue = LineCoordinateArray[CornerPointer];
2287                     LineURValue = LineCoordinateArray[CornerPointer+1];
2288                     LineLLValue = LineCoordinateArray[CornerPointer+ValX];
2289                     if((CornerPointer+1+ValX)<NumberOfCornersInput) {
2290                       LineLRValue = LineCoordinateArray[CornerPointer+1+ValX];
2291                     } else {
2292                       LineLRValue = LineCoordinateArray[CornerPointer+ValX];
2293                     }
2294                     LineMinimum = Math.min(LineULValue,LineURValue);
2295                     LineMinimum = Math.min(LineMinimum,LineLLValue);
2296                     LineMinimum = Math.min(LineMinimum,LineLRValue);
2297    
2298                     /* Test for the presence of a limb in the spline box */
2299    
2300                     if( LineMinimum==-99.9) {
2301                        LABEL30 = true;
2302                     }
2303    
2304                     LineMinimumFinal = (int)(LineMinimum + 0.5);
2305                     if(LineMinimumFinal>MaximumSourceLine) {
2306                        LABEL30 = true;
2307                     }
2308                     LineMaximum = Math.max(LineULValue,LineURValue);
2309                     LineMaximum = Math.max(LineMaximum,LineLLValue);
2310                     LineMaximum = Math.max(LineMaximum,LineLRValue);
2311                     LineMaximumFinal = (int)(LineMaximum + 0.5);
2312                     if(LineMaximumFinal<SourceBlockInc) {
2313                        LABEL30 = true;
2314                     }
2315    
2316                     /* Get 4 corners in elem space & check for out of bound */
2317    
2318                     ValX = remap_vars.ncl;
2319                     ElementULValue = ElementCoordinateArray[CornerPointer];
2320                     ElementURValue = ElementCoordinateArray[CornerPointer+1];
2321                     ElementLLValue = ElementCoordinateArray[CornerPointer+ValX];
2322                     if((CornerPointer+1+ValX)<NumberOfCornersInput) {
2323                       ElementLRValue = ElementCoordinateArray[CornerPointer+1+ValX];
2324                     } else {
2325                       ElementLRValue = ElementCoordinateArray[CornerPointer+ValX];
2326                     }
2327    
2328                     ElementMaximum = Math.max(ElementULValue,ElementURValue);
2329                     ElementMaximum = Math.max(ElementMaximum,ElementLLValue);
2330                     ElementMaximum = Math.max(ElementMaximum,ElementLRValue);
2331                     ElementMaximumFinal = (int)(ElementMaximum + 0.5);
2332                     if(ElementMaximumFinal<1) {
2333                        LABEL30 = true;
2334                     }
2335    
2336                     ElementMinimum = Math.min(ElementULValue,ElementURValue);
2337                     ElementMinimum = Math.min(ElementMinimum,ElementLLValue);
2338                     ElementMinimum = Math.min(ElementMinimum,ElementLRValue);
2339                     ElementMinimumFinal = (int)(ElementMinimum + 0.5);
2340    
2341                     if(ElementMinimumFinal>tiff_vars.in_elems) {
2342                        LABEL30 = true;
2343                     }
2344                     EdgeFixID = 0;
2345    
2346                     /* If the max & min element fall off the image...pitch it */
2347    
2348                     if((ElementMaximumFinal>tiff_vars.in_elems)&&(ElementMinimumFinal<1)) {
2349                        LABEL30 = true;
2350                     }
2351    
2352                     /* Fix if left & right edge should be continuous */
2353    
2354                     if(!LABEL30) { 
2355                        if((ElementMaximumFinal-ElementMinimumFinal)>(int)(.75*tiff_vars.in_elems)) {
2356                           if(ElementULValue<(tiff_vars.in_elems/2)) {
2357                              ElementULValue = ElementULValue+tiff_vars.in_elems;
2358                           }
2359                           if(ElementURValue<(tiff_vars.in_elems/2)) {
2360                              ElementURValue = ElementURValue+tiff_vars.in_elems;
2361                           }
2362                           if(ElementLLValue<(tiff_vars.in_elems/2)) {
2363                              ElementLLValue = ElementLLValue+tiff_vars.in_elems;
2364                           }
2365                           if(ElementLRValue<(tiff_vars.in_elems/2)) {
2366                              ElementLRValue = ElementLRValue+tiff_vars.in_elems;
2367                           }
2368                           EdgeFixID = 1;
2369                        }
2370   
2371                        LineValueA = (LineURValue-LineULValue)/ElementSplineInput;
2372                        LineValueB = (LineLLValue-LineULValue)/LineSplineInput;
2373                        LineValueC = (LineLRValue+LineULValue-LineURValue-LineLLValue)/SplinePixelValue;
2374                        ElementValueA = (ElementURValue-ElementULValue)/ElementSplineInput;
2375                        ElementValueB = (ElementLLValue-ElementULValue)/LineSplineInput;
2376                        ElementValueC = (ElementLRValue+ElementULValue-ElementLLValue-ElementURValue)/SplinePixelValue;
2377   
2378                        int ElementMiscValue = 0;
2379                        LineValueBCounter = LineULValue+0.5;
2380                        LineValueCCounter = 0.0;
2381                        ElementValueBCounter = ElementULValue+0.5;
2382                        ElementValueCCounter = 0.0;
2383   
2384                        if(ReadAllTF==false) {
2385                           ReadAllTF = true;
2386                        }
2387
2388                        if((SplineInc==remap_vars.nspl)||(ElementMinimumFinal<1||(ElementMaximumFinal>tiff_vars.in_elems))||
2389                          ((LineMinimumFinal<SourceBlockInc)||(LineMaximumFinal>MaximumSourceLine))||
2390                          ((FirstCornerPointer+(2*remap_vars.ncl)-1)>LastCornerPointer)) {
2391                           for(YInc=1;YInc<LineSplineInput+1;YInc++) {
2392                              LineValueACounter = LineValueCCounter+LineValueA;
2393                              ElementValueACounter = ElementValueCCounter+ElementValueA;
2394                              LineValueACounter0 = 0.0;
2395                              ElementValueACounter0 = 0.0;
2396                              for(XInc=1;XInc<ElementSplineInput+1;XInc++) {
2397                                 LABEL38 = false;
2398                                 SourceLine = (int)(LineValueBCounter+LineValueACounter0);
2399                                 if(SourceLine<SourceBlockInc) {
2400                                    LABEL38 = true;
2401                                 }
2402                                 if(SourceLine>MaximumSourceLine) {
2403                                    LABEL38 = true;
2404                                 }
2405                                 if(!LABEL38) SourceElement = (int)(ElementValueBCounter+ElementValueACounter0);
2406                                 if(SourceElement<1) {
2407                                    LABEL38 = true;
2408                                 }
2409                                 if((SourceElement>tiff_vars.in_elems)&&(EdgeFixID==0)) {
2410                                    LABEL38 = true;
2411                                 }
2412                                 if(!LABEL38) {
2413                                    if(SourceElement>tiff_vars.in_elems) {
2414                                       SourceElement = SourceElement-tiff_vars.in_elems;
2415                                    }
2416                                    SourceOffsetValue = ((SourceLine-SourceBlockInc)*tiff_vars.in_elems)+SourceElement;
2417                                    OffsetValueE = OffsetValue1+XInc;
2418                                    if(OffsetValueE>tiff_vars.out_elems) {
2419                                       LABEL38 = true;
2420                                    }
2421                                    if(!LABEL38) OffsetValueL = OffsetValue0+ElementMiscValue;
2422                                    if((OffsetValueL/(tiff_vars.out_elems+AccumulatedLines-1))>tiff_vars.out_lines) {
2423                                       LABEL38 = true;
2424                                    }
2425                                    if(!LABEL38) {
2426                                       OffsetValue = OffsetValueL+OffsetValueE;
2427                                       InterpArray[OffsetValue-1] = SourceArray[SourceOffsetValue-1];
2428                                    } /* LABEL38: */
2429                                 } /* LABEL38: */
2430                                 LineValueACounter0 = LineValueACounter0+LineValueACounter;
2431                                 ElementValueACounter0 = ElementValueACounter0+ElementValueACounter;
2432                              }
2433                              LineValueBCounter = LineValueBCounter+LineValueB;
2434                              LineValueCCounter = LineValueCCounter+LineValueC;
2435                              ElementValueBCounter = ElementValueBCounter+ElementValueB;
2436                              ElementValueCCounter = ElementValueCCounter+ElementValueC;
2437                              ElementMiscValue = ElementMiscValue+tiff_vars.out_elems;
2438                           }
2439                        } else {
2440                           if(EdgeFixID==0) {
2441                              for(YInc=1;YInc<LineSplineInput+1;YInc++) {
2442                                  LineValueACounter = LineValueCCounter+LineValueA;
2443                                  ElementValueACounter = ElementValueCCounter+ElementValueA;
2444                                  LineValueACounter0 = 0.0;
2445                                  ElementValueACounter0 = 0.0;
2446                                  OffsetValue = OffsetValue2+ElementMiscValue;
2447                                  for(XInc=1;XInc<ElementSplineInput+1;XInc++) {
2448                                     SourceLine = (int)(LineValueBCounter+LineValueACounter0);
2449                                     SourceElement = (int)(ElementValueBCounter+ElementValueACounter0);
2450                                     SourceOffsetValue = ((SourceLine-SourceBlockInc)*tiff_vars.in_elems)+SourceElement;
2451                                     InterpArray[OffsetValue-1] = SourceArray[SourceOffsetValue-1];
2452                                     OffsetValue = OffsetValue+1;
2453                                     LineValueACounter0 = LineValueACounter0+LineValueACounter;
2454                                     ElementValueACounter0 = ElementValueACounter0+ElementValueACounter;
2455                                  }
2456                                  LineValueBCounter = LineValueBCounter+LineValueB;
2457                                  LineValueCCounter = LineValueCCounter+LineValueC;
2458                                  ElementValueBCounter = ElementValueBCounter+ElementValueB;
2459                                  ElementValueCCounter = ElementValueCCounter+ElementValueC;
2460                                  ElementMiscValue = ElementMiscValue+tiff_vars.out_elems;
2461                              }
2462                           } else if(EdgeFixID==1) {
2463                              for (YInc=1;YInc<LineSplineInput+1;YInc++) {
2464                                 LineValueACounter = LineValueCCounter+LineValueA;
2465                                 ElementValueACounter = ElementValueCCounter+ElementValueA;
2466                                 LineValueACounter0 = 0.0;
2467                                 ElementValueACounter0 = 0.0;
2468                                 OffsetValue = OffsetValue2+ElementMiscValue;
2469                                 for(XInc=1;XInc<ElementSplineInput+1;XInc++) {
2470                                    SourceLine = (int)(LineValueBCounter+LineValueACounter0);
2471                                    SourceElement = (int)(ElementValueBCounter+ElementValueACounter0);
2472                                    if(SourceElement>tiff_vars.in_elems) {
2473                                       SourceElement = SourceElement-tiff_vars.in_elems;
2474                                    }
2475                                    SourceOffsetValue = ((SourceLine-SourceBlockInc)*tiff_vars.in_elems)+SourceElement;
2476                                    InterpArray[OffsetValue-1] = SourceArray[SourceOffsetValue-1];
2477                                    OffsetValue = OffsetValue+1;
2478                                    LineValueACounter0 = LineValueACounter0+LineValueACounter;
2479                                    ElementValueACounter0 = ElementValueACounter0+ElementValueACounter;
2480                                 }
2481                                 LineValueBCounter = LineValueBCounter+LineValueB;
2482                                 LineValueCCounter = LineValueCCounter+LineValueC;
2483                                 ElementValueBCounter = ElementValueBCounter+ElementValueB;
2484                                 ElementValueCCounter = ElementValueCCounter+ElementValueC;
2485                                 ElementMiscValue = ElementMiscValue+tiff_vars.out_elems;
2486                              }
2487                           }
2488                        }
2489                     } /* LABEL30: */
2490                     SourceBlockValue = SourceBlockValue+1;
2491                  }
2492                  FirstCornerPointer = FirstCornerPointer+remap_vars.ncl;
2493               }
2494            }
2495            BlockCounter = BlockCounter+1;
2496         }
2497             
2498         for(YInc=0;YInc<tiff_vars.out_lines;YInc++) {
2499            for(XInc=0;XInc<tiff_vars.out_elems;XInc++) {
2500               ArrayIndexValue = (YInc*tiff_vars.out_elems)+XInc;
2501               if((ArrayIndexValue>0)&&(InterpArray[ArrayIndexValue]==0.0)) {
2502                  InterpArray[ArrayIndexValue] = InterpArray[ArrayIndexValue-1];
2503               }
2504               IRData_Remap_Temperature[YInc][XInc] = InterpArray[ArrayIndexValue];
2505            }
2506         }
2507         Ret_ID = 0;
2508      }
2509
2510      return Ret_ID;
2511
2512   }
2513
2514   private static int IND(int y, int x) {
2515       /* #define IND(y,x) ((y-1)*remap_vars.ncl+x) */
2516      return ((y-1)*remap_vars.ncl+x);
2517   } 
2518
2519}