/*
VisAD system for interactive analysis and visualization of numerical
data.  Copyright (C) 1996 - 2017 Bill Hibbard, Curtis Rueden, Tom
Rink, Dave Glowacki, Steve Emmerson, Tom Whittaker, Don Murray, and
Tommy Jasmin.

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Library General Public License for more details.

You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA
*/

/*
RH looks interesting up to 10000 m
MR better near boundary layer (but tighter color range)
*/

package visad.aeri;

import visad.*;
import visad.util.*;
import visad.DisplayImpl;
import visad.java3d.DisplayImplJ3D;
import visad.data.netcdf.*;
import visad.bom.WindPolarCoordinateSystem;
import visad.data.mcidas.BaseMapAdapter;
import visad.data.mcidas.AreaAdapter;
import visad.meteorology.ImageSequenceManager;
import visad.meteorology.NavigatedImage;
import visad.meteorology.ImageSequence;
import java.rmi.RemoteException;
import java.io.IOException;
import java.io.File;
import visad.data.visad.VisADSerialForm;

// JFC packages
import javax.swing.*;

// AWT packages
import java.awt.*;
import java.awt.event.*;

public class Aeri
       implements ScalarMapListener
{
  RealType latitude;
  RealType longitude;
  RealType altitude;

  //- (lon,lat,alt)
  //
  RealTupleType spatial_domain;
  RealType time;
  RealType stn_idx;

  RealType temp;
  RealType dwpt;
  RealType wvmr;
  RealType RH;
  RealType theta;
  RealType thetaE;
  RealType band1 = null;

  //- (T,TD,WV,AGE)
  //
  RealTupleType advect_range;

  FunctionType advect_type;
  FunctionType advect_field_type;

  FieldImpl advect_field;
  FieldImpl stations_field;
  ImageSequence image_seq;

  int n_stations = 5;

  double[] station_lat;
  double[] station_lon;
  double[] station_alt;
  double[] station_id;

  BaseMapAdapter baseMap;
  DataReference map_ref;
  DataReference img_ref;

  ScalarMap xmap;
  ScalarMap ymap;
  ScalarMap zmap;
  ScalarMap img_map;
  boolean xmapEvent = false;
  boolean ymapEvent = false;
  boolean imgEvent = false;
  boolean first_xy_Event = false;
  boolean first_img_Event = false;

  float latmin, latmax;
  float lonmin, lonmax;
  float del_lat, del_lon;
  double[] x_range, y_range;

  int height_limit = 3000; // meters
  boolean rh = false;
  boolean tm = false;
  boolean pt = false;
  boolean ept = false;

  int start_date = 0;
  double start_time = 0.0;

  public static void main(String args[])
         throws VisADException, RemoteException, IOException
  {
    Aeri aeri = new Aeri(args);
  }

  public Aeri(String[] args)
         throws VisADException, RemoteException, IOException
  {
    String vadfile = null;
    String baseDate = "20000112";
    for (int i=0; i hgt) hgt = hi[2];
        if (!any && samples[0][0] == samples[0][0] &&
            samples[1][0] == samples[1][0]) {
          any = true;

          locs[0][0] = samples[0][0];
          locs[1][0] = samples[1][0];
          locs[0][1] = samples[0][0];
          locs[1][1] = samples[1][0];

          if (samples[0][0] > lat_max) lat_max = samples[0][0];
          if (samples[0][0] < lat_min) lat_min = samples[0][0];
          if (samples[1][0] > lon_max) lon_max = samples[1][0];
          if (samples[1][0] < lon_min) lon_min = samples[1][0];
        }
      }
      if (any) {
        locs[2][0] = 0.0f;
        locs[2][1] = hgt;
        set_s[ii++] = new Gridded3DSet(spatial_domain, locs, 2, null, null, null);
// System.out.println("set_s[" + kk + "] = " + set_s[kk]);
        if (hgt > hgt_max) hgt_max = hgt;
      }
    }
    SampledSet[] set_ss = new SampledSet[ii];
    System.arraycopy(set_s, 0, set_ss, 0, ii);
    return new UnionSet(spatial_domain, set_ss);
  }

  public void mapChanged(ScalarMapEvent e)
       throws VisADException, RemoteException
  {
    if ( xmap.equals(e.getScalarMap()) ) {
      xmapEvent = true;
    }
    else if ( ymap.equals(e.getScalarMap()) ) {
      ymapEvent = true;
    }
    else if ( img_map.equals(e.getScalarMap()) ) {
      imgEvent = true;
    }
    if (( xmapEvent && ymapEvent ) && !(first_xy_Event) ) {
      x_range = xmap.getRange();
      y_range = ymap.getRange();
      latmin = (float)y_range[0];
      latmax = (float)y_range[1];
      lonmin = (float)x_range[0];
      lonmax = (float)x_range[1];
      baseMap.setLatLonLimits(latmin-del_lat, latmax+del_lat,
                              lonmin-del_lon, lonmax+del_lon);
      DataImpl map = baseMap.getData();
      map_ref.setData(map);
      first_xy_Event = true;
      xmap.setRange(lonmin, lonmax);
      ymap.setRange(latmin, latmax);
      img_ref.setData(image_seq);
    }
    if ( imgEvent && !first_img_Event ) {
      double[] i_range = img_map.getRange();
      System.out.println(i_range[0]+" "+i_range[1]);
      first_img_Event = true;
      img_map.setRange(i_range[1], i_range[0]);
    }
  }

  public void controlChanged(ScalarMapControlEvent e)
  {
  }

  FieldImpl[] makeWinds(String[] files)
              throws VisADException, RemoteException, IOException
  {
    DataImpl[] file_data = new DataImpl[n_stations];
    FieldImpl[] time_field = new FieldImpl[n_stations];
    double[][] time_offset = null;
    double[] base_date = new double[n_stations];
    double[] base_time = new double[n_stations];
    Gridded1DSet d_set = null;
    FlatField new_ff = null;

    RealType alt;
    RealType spd;
    RealType dir;

    RealType u_wind = RealType.getRealType("u_wind");
    RealType v_wind = RealType.getRealType("v_wind");

    //- create a new netcdf reader
    Plain plain = new Plain();

    //- retrieve file data objects
    for ( int kk = 0; kk < n_stations; kk++) {
      file_data[kk] = plain.open(files[kk]);
    }

    //- make sub mathtype for file objects
    MathType file_type = file_data[0].getType();
    System.out.println(file_type.prettyString());
    System.out.println();

    FunctionType f_type0 =
      (FunctionType)((TupleType)file_type).getComponent(2);

    TupleType ttype = (TupleType) f_type0.getRange();
    int p_field_index = ttype.getDimension() - 1;
    FunctionType f_type1 =
  //- (FunctionType)((TupleType)f_type0.getRange()).getComponent(10);
  //- (FunctionType)((TupleType)f_type0.getRange()).getComponent(12);
      (FunctionType) ttype.getComponent(p_field_index);

    RealTupleType rng_tuple = (RealTupleType) f_type1.getRange();
    int altitude_index = rng_tuple.getIndex("Altitude");
 //-altitude = (RealType)((RealTupleType)f_type1.getRange()).getComponent(0);
 // altitude = (RealType)((RealTupleType)f_type1.getRange()).getComponent(1);
    altitude = (RealType) rng_tuple.getComponent(altitude_index);
    int windSpeed_index = rng_tuple.getIndex("windSpeed");
 //-spd = (RealType)((RealTupleType)f_type1.getRange()).getComponent(4);
 //-spd = (RealType)((RealTupleType)f_type1.getRange()).getComponent(5);
    spd = (RealType) rng_tuple.getComponent(windSpeed_index);
    int windDir_index = rng_tuple.getIndex("windDir");
 //-dir = (RealType)((RealTupleType)f_type1.getRange()).getComponent(3);
 //-dir = (RealType)((RealTupleType)f_type1.getRange()).getComponent(4);
    dir = (RealType) rng_tuple.getComponent(windDir_index);
    RealType[] r_types = { dir, spd };

    /* WLH 28 Dec 99 */
    RealType[] uv_types = { u_wind, v_wind };
    // EarthVectorType uv = new EarthVectorType(uv_types); WLH meeds m/s
    RealTupleType uv = new RealTupleType(uv_types);
    CoordinateSystem cs = new WindPolarCoordinateSystem(uv);
    RealTupleType ds = new RealTupleType(r_types, cs, null);


    FunctionType alt_to_ds = new FunctionType(altitude, ds);
    FunctionType alt_to_uv = new FunctionType(altitude, uv);

    RealType domain_type = (RealType)
      ((TupleType)f_type0.getRange()).getComponent(0);
    time = domain_type;
    FunctionType new_type = new FunctionType(RealType.Time, alt_to_uv);

    FieldImpl[] winds = new FieldImpl[n_stations];

    for ( int ii = 0; ii < n_stations; ii++ )
    {
      base_date[ii] = (double)
        ((Real)((Tuple)file_data[ii]).getComponent(0)).getValue();
      base_time[ii] = (double)
        ((Real)((Tuple)file_data[ii]).getComponent(1)).getValue();
      time_field[ii] = (FieldImpl)
        ((Tuple)file_data[ii]).getComponent(2);
      station_lat[ii] =
        ((Real)((Tuple)time_field[ii].getSample(0)).getComponent(6)).getValue();
      station_lon[ii] =
       -((Real)((Tuple)time_field[ii].getSample(0)).getComponent(7)).getValue();
      station_alt[ii] =
        ((Real)((Tuple)time_field[ii].getSample(0)).getComponent(8)).getValue();
      station_id[ii] =
        ((Real)((Tuple)time_field[ii].getSample(0)).getComponent(9)).getValue();
      if (ii == 0) {
        start_time = base_time[0];
        start_date = (int) base_date[0];
      }

/*
System.out.println("wind " + ii + " date: " + base_date[ii] +
                   " time: " + base_time[ii]);
wind 0 date: 1.9991226E7 time: 9.461664E8
wind 1 date: 1.9991226E7 time: 9.461664E8
wind 2 date: 1.9991226E7 time: 9.461664E8
wind 3 date: 1.9991226E7 time: 9.461664E8
wind 4 date: 1.9991226E7 time: 9.461664E8

looks like 19991226 and seconds since 1-1-70
time_offset looks like seconds since 0Z
*/

      int length = time_field[ii].getLength();
      time_offset = new double[1][length];
      FlatField[] range_data = new FlatField[length];

      double[][] samples = null; // WLH
      for ( int jj = 0; jj < length; jj++ )
      {
        Tuple range = (Tuple) time_field[ii].getSample(jj);
        time_offset[0][jj] = (double)((Real)range.getComponent(0)).getValue();

     //-FlatField p_field = (FlatField) range.getComponent(10);
     //-FlatField p_field = (FlatField) range.getComponent(12);
        FlatField p_field = (FlatField) range.getComponent(p_field_index);
        double[][] values =
          p_field.getValues(); // WLH - (alt, ?, ?, dir, spd, ???)
        double[][] new_values = new double[2][values[0].length];

        if ( jj == 0 )  //- only once, vertical range gates don't change
        {
          samples = new double[1][values[0].length]; // WLH
       //-System.arraycopy(values[0], 0, samples[0], 0, samples[0].length);
       //-System.arraycopy(values[1], 0, samples[0], 0, samples[0].length);
          System.arraycopy(values[altitude_index], 0, samples[0], 0, samples[0].length);
          d_set = new Gridded1DSet(altitude, Set.doubleToFloat(samples), samples[0].length);
        }
        new_ff = new FlatField(alt_to_uv, d_set);

/* WLH - fill in missing winds - this also extrapolates to missing
start or end altitudes - but it does nothing if all winds are missing */
        int n_not_miss = 0;
        int[] not_miss = new int[values[0].length];
        for ( int mm = 0; mm < values[0].length; mm++ )
        {
       //-if ( values[3][mm] <= -9999 ) {
       //-if ( values[4][mm] <= -9999 ) {
          if ( values[windDir_index][mm] <= -9999 ) {
            new_values[0][mm] = Float.NaN;
          }
          else {
         //-new_values[0][mm] = values[3][mm];
         //-new_values[0][mm] = values[4][mm];
            new_values[0][mm] = values[windDir_index][mm];
          }

       //-if ( values[4][mm] <= -9999 ) {
       //-if ( values[5][mm] <= -9999 ) {
          if ( values[windSpeed_index][mm] <= -9999 ) {
            new_values[1][mm] = Float.NaN;
          }
          else {
         //-new_values[1][mm] = values[4][mm];
         //-new_values[1][mm] = values[5][mm];
            new_values[1][mm] = values[windSpeed_index][mm];
          }

          if (new_values[0][mm] == new_values[0][mm] &&
              new_values[1][mm] == new_values[1][mm]) {
            not_miss[n_not_miss] = mm;
            n_not_miss++;
          }
        }
        if (0 < n_not_miss && n_not_miss < values[0].length) {
          int nn = n_not_miss;
          if (not_miss[0] > 0) nn += not_miss[0];
          int endlen = values[0].length - (not_miss[n_not_miss-1]+1);
          if (endlen > 0) nn += endlen;

          float[][] newer_values = new float[2][nn];
          float[][] newer_samples = new float[1][nn];
          // fill in non-missing values
          for (int i=0; i (lon,lat,alt) -> (T,TD,WV,AGE)
    //
    advect_field = new FieldImpl(advect_field_type, rtvls.getDomainSet());

    //- resample winds domain (time) to rtvls
    //

    rtvls_domain = rtvls.getDomainSet();
    wind_to_rtvl_time = (FieldImpl)
      winds.resample( rtvls_domain,
                      Data.WEIGHTED_AVERAGE,
                      Data.NO_ERRORS );


    //- resample rtvls domain (altitude) to winds at each rtvl time
    //
    int len = rtvls.getLength();
    rtvl_on_wind = new FlatField[len];
    wind_on_wind = new FieldImpl[len];
    for ( int tt = 0; tt < len; tt++ ) {
      alt_to_rtvl = (FieldImpl) rtvls.getSample(tt);
      alt_to_wind = (FieldImpl) wind_to_rtvl_time.getSample(tt);
      Set ds = alt_to_wind.getDomainSet();
// WLH height limit
      float[][] samples = ds.getSamples();
      float[] ns = new float[samples[0].length];
      int nn = 0;
      for (int i=0; i -50.) {
      retval = ( coef[0] + temp*(coef[1] + temp*(coef[2] + temp*(coef[3] +
      temp*(coef[4] + temp*(coef[5] + temp*(coef[6] + temp*coef[7])))))) );
    }
    else {
       double tt = (-temp - 50.)/5.;
       int inx = (int) tt;
       if (inx < escold.length) {
         retval = escold[inx] + (tt % 1.)*(escold[inx+1]-escold[inx]);
       } else {
         retval = 1e-7;
       }
    }
    return retval;
  }

  /** mixing ratio
  *
  */
  public static double mixingRatio(double t, double p) {
    double e = satVapPres(t);
    return ( 621.97*e/(p - e) );
  }

  /** relative humidity
  *
  */
  public static double relativeHumidity(double t, double td) {
    return (satVapPres(td) / satVapPres(t));
  }

  public static double potentialTemperature( double t, double p ) 
  {
    double R_Cp = 0.28585;
    double Ps = 1000;  //- reference pressure, (mb)

    return t*Math.pow((Ps/p), R_Cp);
  }

  public static double equivPotentialTemperature(double theta, 
                                                 double temp,
                                                 double press )
  {
    double Lc = 2.5*1000000.0;
    double Cp = 1004.0;
    double thetaE;
    double Qs;

    Qs = mixingRatio(temp, press*100);
    Qs = Qs/1000.0;
    thetaE = theta*Math.exp((Lc*Qs)/(Cp*temp));
    return thetaE;
  }

  public static double equivPotentialTemperatureStar( double theta,
                                                      double Q,
                                                      double temp )
  {
    double Lc = 2.5*1000000.0;
    double Cp = 1004.0;

    Q = Q/1000.0;
    return theta*Math.exp((Lc*Q)/(Cp*temp));
  }
}