//////////////////////////////////////////////////
// JIST (Java In Simulation Time) Project
// Timestamp: <Fading.java Wed 2004/06/23 09:16:53 barr pompom.cs.cornell.edu>
//

// Copyright (C) 2004 by Cornell University
// All rights reserved.
// Refer to LICENSE for terms and conditions of use.

package jist.swans.field;

import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;
import java.util.StringTokenizer;

import org.apache.commons.math.stat.StatUtils;
import jist.runtime.JistAPI;
import jist.runtime.Main;
import jist.swans.Constants;
import jist.swans.misc.Location;
import jist.swans.misc.Util;
import jist.swans.radio.RadioInfo;

/**
 * Interface for performing fading calculations.
 *
 * @author Rimon Barr &lt;barr+jist@cs.cornell.edu&gt;
 * @version $Id: Fading.java,v 1.10 2004/06/23 17:15:54 barr Exp $
 * @since SWANS1.0
 */

public interface Fading
{
  //////////////////////////////////////////////////
  // interface
  //

  /**
   * Compute the fading loss.
   *
   * @deprecated
   * @return fading loss (units: dB)
   */
  double compute();

  /**
   * Compute the fading loss.
   *
   * @param srcRadio source radio information
   * @param srcLocation source location
   * @param dstRadio destination radio information
   * @param dstLocation destination location
   * @return fading loss (units: dB)
   */
  double compute(RadioInfo srcRadio, Location srcLocation,
      RadioInfo dstRadio, Location dstLocation);

  /**
   * Compute the max. signal strength between source radio and given point.
   *
   * @param srcRadio source radio information
   * @param srcLocation source location
   * @param dstRadio destination radio information
   * @param dstLocation destination location
   * @return fading loss (units: dB)
   */
  double computeLimit(RadioInfo srcRadio, Location srcLocation,
      RadioInfo dstRadio, Location dstLocation);

  //////////////////////////////////////////////////
  // implementations
  //

  /**
   * Computes zero fading.
   *
   * @author Rimon Barr &lt;barr+jist@cs.cornell.edu&gt;
   * @since SWANS1.0
   */
  final class None implements Fading
  {
    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute()
     */
    public double compute(){
      return .0;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double compute(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      return .0;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double computeLimit(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      return .0;
    }

  }

  /**
   * Computes Rayleigh fading. Equivalent to GloMoSim code.
   *
   * contructs a rayleigh distribution with sigma = 0.7979 = 1/sqrt(pi/2)
   *
   * TODO if X,Y ~ N(0,sigma^2), then R=sqrt(X^2+Y^2) ~ Rayleigh(sigma) and
   * R^2 ~ Exponential, where 2 * sigma^2 = Power(pathloss+shadowing)
   *
   * The Rayleigh distribution has the cdf as:
   * F = 1 - exp(-r^2/(2*sigma^2))
   *
   * So we generate r using inverse transform method as:
   * r = sqrt(-2*sigma^2*log(u)),
   * where u is a uniform random number between 0 and 1.
   *
   * To compute the channel gain, we need to compute:
   * fading_gain = 10*log10(r^2) = 10*log10(-2*sigma^2*log(u)).
   * However, the original code is 5*log10(-2*sigma^2*log(u)).
   * For the Rayleigh variance, the expected value of r^2 is:
   * E(r^2) = 2*sigma^2, we want this value to be 1, thus sigma^2 has to be 1/2.
   *
   * @author Rimon Barr &lt;barr+jist@cs.cornell.edu&gt;
   * @since SWANS1.0
   */
  final class Rayleigh implements Fading
  {
    /** Rayleigh distribution variance constant. */
    private static final double VARIANCE = 0.6366197723676; // = 1/sqrt(pi/2)

    // Fading interface
    /** {@inheritDoc} */
    public double compute()
    {
      // compute fading_dB; positive values are signal gains
      return 5.0 * Math.log(-2.0 * VARIANCE * Math.log(Constants.random.nextDouble())) / Constants.log10;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double compute(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      return compute();
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double computeLimit(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      throw new RuntimeException("not implemented");
    }
  }

  /**
   * Computes Rician fading. Equivalent to GloMoSim code.
   *
   * @author Rimon Barr &lt;barr+jist@cs.cornell.edu&gt;
   * @since SWANS1.0
   */
  final class Rician implements Fading
  {
    /** distribution parameters. */
    private final double kFactor, stddev;

    /**
     * Create new Rician fading model object.
     *
     * @param kFactor k
     */
    public Rician(double kFactor)
    {
      this.kFactor = kFactor;
      this.stddev = computeStandardDeviation(kFactor);
    }

    /**
     * Compute zero-order Bessel function.
     *
     * @param x input
     * @return output of Bessel
     */
    private static double Besseli0(double x)
    {
      double ax = Math.abs(x);
      if (ax < 3.75)
      {
        double y = x/3.75;
        y *= y;
        return 1.0 + y*(3.5156229 + y*(3.0899424 + y*(1.2067492 +
          y*(0.2659732 + y*(0.360768e-1 + y*0.45813e-2)))));
      }
      else
      {
        double y = 3.75/ax;
        return (Math.exp(ax)/Math.sqrt(ax)) * (0.39894228 + y*(0.1328592e-1 +
          y*(0.225319e-2 + y*(-0.157565e-2 + y*(0.916281e-2 + y*(-0.2057706e-1 +
          y*(0.2635537e-1 + y*(-0.1647633e-1 + y*0.392377e-2))))))));
      }
    }

    /**
     * Compute first-order Bessel function.
     *
     * @param x input
     * @return output of Bessel
     */
    private static double Besseli1(double x)
    {
      double ax = Math.abs(x);
      if (ax < 3.75)
      {
        double y = x/3.75;
        y *= y;
        return x * (0.5 + y*(0.87890494 + y*(0.51498869 + y*(0.15084934 + y*(0.2658733e-1 +
          y*(0.301532e-2 + y*0.32411e-3))))));
      }
      else
      {
        double y = 3.75/ax;
        return Math.abs((Math.exp(ax)/Math.sqrt(ax)) * (0.39894228 + y*(-0.3988024e-1 +
          y*(-0.362018e-2 + y*(0.163801e-2 + y*(-0.1031555e-1 + y*(0.2282967e-1 +
          y*(-0.2895312e-1 + y*(0.1787654e-1 - y*0.420059e-2)))))))));
      }
    }

    /**
     * Computes standard deviation for Rician distribution such that mean is 1.
     *
     * @param kFactor k
     * @return Rician standard deviation
     */
    private static double computeStandardDeviation(double kFactor)
    {
      return 1.0/(Math.sqrt(Math.PI/2.0) * Math.exp(-kFactor/2.0) *
        ((1+kFactor)*Besseli0(kFactor/2.0) + kFactor*Besseli1(kFactor/2.0)));
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute()
     */
    public double compute()
    {
      // compute fading_dB; positive values are signal gains
      double a = Math.sqrt(2.0 * kFactor * stddev * stddev), r, v1, v2;
      do
      {
        v1 = -1.0 + 2.0 * Constants.random.nextDouble();
        v2 = -1.0 + 2.0 * Constants.random.nextDouble();
        r = v1 * v1 + v2 * v2;
      }
      while (r >= 1.0 || r == 0.0); // see java.util.Random.nextGaussian()
      r = Math.sqrt(-2.0 * Math.log(r) / r);
      v1 = a + stddev * v1 * r;
      v2 = stddev * v2 * r;
      return 5.0 * Math.log(v1*v1 + v2*v2) / Constants.log10;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double compute(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      return compute();
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double computeLimit(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      throw new RuntimeException("not implemented");
    }
  }


  /**
   * Precalculated Rician fading according to
   * Punnoose "Efficient Simulation of Ricean Fading within a Packet Simulator"
   * revised according to http://www.winlab.rutgers.edu/~zhibinwu/html/ns_fading_error.html
   *
   * TODO should we use randomized time offsets per node pair??
   *
   * maxVelocity
   *   This parameter determines the doppler spectrum.
   *   The value of this parameter is the maximum speed of any objects in the
   *   environment (could be nodes).  It is NOT the speed of the
   *   nodes themselves.
   *
   * K
   *   Ricean K factor
   *
   *
   * Doppler spreas fm = maxVelocity / lambda = maxVelocity / .12 [2.4GHz]
   *
   * Coherence Time Tc = k / fm,  k = .25 ... .5
   *   e.g. for maxVelocity=.1  -> Tc=300ms
   *   e.g. for maxVelocity=1   -> Tc=30ms
   *   e.g. for maxVelocity=2.5 -> Tc=12ms
   *   e.g. for maxVelocity=5   -> Tc=6ms
   *   e.g. for maxVelocity=10  -> Tc=3ms
   *
   * For outdoor LOS environments, do not use fading. Use correlated shadowing
   * with shadowing coherence 1-10 sec (see
   * http://fly.isti.cnr.it/curriculum/papers/pdf/Rural-model-Winmee07.pdf).
   *
   * @author kurth
   */
  public static final class PunnooseRician implements Fading {

    /** Maximum velocity of vehicle/objects in environment.
     * Used for computing doppler */
    protected double maxVelocity = 2.5;

    /** Ricean K factor */
    protected double K = Util.fromDB(6);

    /** Data values for inphase and quad phase */
    protected double[][] data = new double[2][];
    /** Max doppler freq in table */
    private double fm0;
    /** Sampling rate */
    protected double fs;

    protected int[][] startIdx;

    private boolean symmetric = true;

    /** max fading power at the given percentile, calculated from data file */
    private double fadeMargin;

    private double percentile;


    public PunnooseRician(double maxVelocity, double riceanKdB, //double lambda,
        String traceFile, boolean symmetric) throws IOException {
      this.maxVelocity = maxVelocity;
      this.K = Util.fromDB(riceanKdB);
      this.symmetric = symmetric;
      this.percentile = 99.;
      if (Main.ASSERT)
        Util.assertion(!Double.isInfinite(this.K) && !Double.isNaN(this.K));
      loadDataFile(traceFile);

      int N = data[0].length;
      startIdx = new int[20][20];
      for (int i = 0; i < startIdx.length; i++)
        for (int j = 0; j < startIdx[i].length; j++)
          startIdx[i][j] = Constants.random.nextInt(N);
    }

    private void loadDataFile(String traceFile) throws IOException {
      BufferedReader reader = new BufferedReader(new FileReader(traceFile));
      int N = -1;

      while(true) {
        String line = reader.readLine();
        if (line.trim().startsWith("#"))
          continue;

        StringTokenizer tok = new StringTokenizer(line);
        String arg1 = tok.nextToken();

        if (arg1.equals("DATA"))
          break;

        if(arg1.equals("fm")) {
          fm0 = Double.parseDouble(tok.nextToken());
        } else if(arg1.equals("fs")) {
          fs = Double.parseDouble(tok.nextToken());
        } else if(arg1.equals("N")) {
          N = Integer.parseInt(tok.nextToken());
        }
      }

      data[0] = new double[N];
      data[1] = new double[N];

      double[] power_dB = new double[N];
      for(int k=0; k<N; k++) {
        StringTokenizer tok = new StringTokenizer(reader.readLine());
        data[0][k] = Double.parseDouble(tok.nextToken());
        data[1][k] = Double.parseDouble(tok.nextToken());

        // Find the envelope multiplicative factor
        double tmp_x1 = data[0][k] + Math.sqrt(2.0 * K);
        double tmp_x2 = data[1][k];

        power_dB[k] = Util.toDB((tmp_x1*tmp_x1 + tmp_x2*tmp_x2) / (2.0 * (K+1)));
      }

      this.fadeMargin = StatUtils.percentile(power_dB, percentile);
//      for (double i = .01; i < 1.; i += 0.01)
//        System.out.println(i + "       " + StatUtils.percentile(power_dB, i*100));
    }

    /**
     *
     * @return fading power factor (absolute, not dB!)
     */
    private final double Pr(double lambda, int idxTr, int idxRx) {
      // Compute a unique table offset for unique pairs of nodes
      int N = data[0].length;

      // Max doppler freq in scenario
      double fm = maxVelocity / lambda;

      double timeNow = JistAPI.getTime() / (double)Constants.SECOND;
      double time_index = (timeNow) * fs * fm / fm0;

      if (this.symmetric ) {
        int max = Math.max(idxTr, idxRx);
        int min = Math.min(idxTr, idxRx);
        idxTr = max;
        idxRx = min;
      }

      try {
        // use randomized time offsets per node pair
        time_index = (time_index + startIdx[idxTr][idxRx]) % N;
      } catch (ArrayIndexOutOfBoundsException e) {
        int oldLen = startIdx.length;
        int newLen = Math.max(idxTr, idxRx) + 17;
        int[][] newIdx = new int[newLen][];
        System.arraycopy(startIdx, 0, newIdx, 0, oldLen);
        startIdx = newIdx;
        for (int i = 0; i < startIdx.length; i++) {
          int start = 0;
          int[] row = new int[newLen];
          if (i < oldLen) {
            System.arraycopy(startIdx[i], 0, row, 0, oldLen);
            start = oldLen;
          }
          startIdx[i] = row;
          for (int j = start; j < startIdx[i].length; j++)
            startIdx[i][j] = Constants.random.nextInt(N);
        }
        time_index = (time_index + startIdx[idxTr][idxRx]) % N;
      }

      // Do envelope interpolation using Legendre polynomials
      double X0, X1, X2, X3;
      int ind0, ind1, ind2, ind3;

      ind1 = (int) Math.floor(time_index);
      ind0 = (ind1-1+N) % N;
      ind2 = (ind1+1) % N;
      ind3 = (ind1+2) % N;

      X1 = time_index - ind1;
      X0 = X1+1.0;
      X2 = X1-1.0;
      X3 = X1-2.0;

      double x1_interp = data[0][ind0]*X1*X2*X3/(-6.0) +
        data[0][ind1]*X0*X2*X3*(0.5) +
        data[0][ind2]*X0*X1*X3*(-0.5) +
        data[0][ind3]*X0*X1*X2/6.0;

      double x2_interp = data[1][ind0]*X1*X2*X3/(-6.0) +
        data[1][ind1]*X0*X2*X3*(0.5) +
        data[1][ind2]*X0*X1*X3*(-0.5) +
        data[1][ind3]*X0*X1*X2/6.0;

      // Find the envelope multiplicative factor
      double tmp_x1 = x1_interp + Math.sqrt(2.0 * K);
      double tmp_x2 = x2_interp;

      double envelope_fac = (tmp_x1*tmp_x1 + tmp_x2*tmp_x2) / (2.0 * (K+1));

      //Pr_Rice = 10.0 * log10(envelope_fac);
      double Pr_Rice = envelope_fac;

      return Pr_Rice;
    }


    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute()
     */
    public double compute() {
      // needs location information
      throw new UnsupportedOperationException();
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double compute(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      return Util.toDB(Pr(srcRadio.getWaveLength(), srcRadio.getId(), dstRadio.getId()));
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double computeLimit(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      return fadeMargin;
    }
  }


  /**
   * Propagation Model from http://www.dei.unipd.it/wdyn/?IDsezione=5529
   *
   * @author Federico Maguolo <maguolof@dei.unipd.it>
   * @author kurth
   */
  public static class Jakes implements Fading {
    protected static class ComplexNumber {
      public double real;
      public double imag;

      public ComplexNumber(double real, double imag) {
        super();
        this.real = real;
        this.imag = imag;
      }
    }

    protected class PathCoefficients {
      private RadioInfo receiver;
      private int nOscillators;
      private int nRays;
      private double[][] phases;
      private long start;

      public PathCoefficients (RadioInfo receiver, int nRays, int nOscillators) {
        this.receiver = receiver;
        this.nRays = nRays;
        this.nOscillators = nOscillators;
        DoConstruct ();
      }
      private void DoConstruct () {
        phases = new double[nRays][];
        for (int i = 0; i < nRays; i++) {
          phases[i] = new double[nOscillators + 1];
          for (int j = 0; j <= nOscillators; j++) {
            phases[i][j] = 2.0 * Math.PI * Jakes.this.variable;
          }
        }
        start = JistAPI.getTime();
      }
      public RadioInfo getReceiver () {
        return receiver;
      }

      /**
       *
       * @param txPowerDbm
       * @return rxPower or fading coefficient [in dB]
       */
      public double getRxPowerDbm (double txPowerDbm) {
        int N = 4 * nOscillators + 2;
        double interval = JistAPI.getTime() - start;
        ComplexNumber fading = new ComplexNumber(0.0, 0.0);
        for (int i = 0; i < nRays; i++) {
          for (int j = 0; j <= nOscillators; j++) {
            double phase = phases[i][j] + 2.0 * Math.PI * Math.cos (2.0 * Math.PI * j / N)
              * Jakes.this.fd * (interval / Constants.SECOND);
            fading.real += Jakes.this.amp[j].real * Math.cos (phase);
            fading.imag += Jakes.this.amp[j].imag * Math.cos (phase);
          }
        }
        double coef = Math.sqrt(Math.pow (fading.real, 2.0) + Math.pow (fading.imag, 2.0));
        return txPowerDbm + Util.toDB(coef);
      }
    }

    private static class PathsSet {
      public RadioInfo sender;
      public List<PathCoefficients> receivers;
      public PathsSet(RadioInfo sender) {
        this.sender = sender;
        this.receivers = new ArrayList<PathCoefficients>();
      }
    }

    private List<PathsSet> paths;

    protected ComplexNumber[] amp;

    /**
     * JakesPropagationNumberOfRaysPerPath - The number of rays to use for
     * compute the fading coeficent for a given path (default is 1)
     */
    private int nRays = 1;

    /**
     * JakesPropagationNumberOfOscillatorsPerRay - The number of oscillators to
     * use for compute the coeficent for a given ray of a given path (default is 4)
     */
    private int nOscillators = 4;


//    /**
//     * JakesPropagationLossDistribution - The distribution to choose the initial
//     * phases (default Constant:1.0)
//     */
//    protected Random variable;
    protected int variable = 1;

    /**
     * JakesPropagationDopplerFreq - The doppler frequency in Hz
     * (f_d = v / lambda = v * f / c, the defualt is 0)
     */
    protected double fd = 0;

    public Jakes () {
      DoConstruct();
    }

    public Jakes(int rays, int oscillators, double fd) {
      nRays = rays;
      nOscillators = oscillators;
      this.fd = fd;
      DoConstruct();
    }

//    public Jakes (RandomVariable variable) {
//
//    }

    private void DoConstruct () {
      int N = 4 * nOscillators + 2;
      amp = new ComplexNumber[nOscillators + 1];
      amp[0] = new ComplexNumber(2.0 * Math.sqrt(2.0 / N) * Math.cos (Math.PI / 4.0),
          2.0 * Math.sqrt(2.0 / N) * Math.sin (Math.PI / 4.0));
      for (int i = 1; i <= nOscillators; i++) {
        double beta = Math.PI * (double)i / nOscillators;
        amp[i] = new ComplexNumber(4.0 * Math.cos (beta) / Math.sqrt(N),
            4.0 * Math.sin (beta) / Math.sqrt(N));
      }
      this.paths = new ArrayList<PathsSet>();
    }

    public void setNRays (int nRays) {
      this.nRays = nRays;
    }

    public void setNOscillators (int nOscillators) {
      this.nOscillators = nOscillators;
    }

    /**
    *
    * @param txPowerDbm
    * @return rxPower or fading coefficient [in dB]
    */
    public double getRxPower (double txPowerDbm, RadioInfo a, RadioInfo b) {
      for (int i = paths.size()-1; i >= 0; i--) {
        PathsSet ps = paths.get(i);
        if (ps.sender.equals(a)) {
          paths.remove(i);
          paths.add(ps);
          Iterator iter = ps.receivers.iterator();
          while (iter != null && iter.hasNext()) {
            PathCoefficients pc = (PathCoefficients) iter.next();
            if (pc.getReceiver().equals(b)) {
              ps.receivers.remove(iter);
              ps.receivers.add(pc);
              return pc.getRxPowerDbm (txPowerDbm);
            }
          }
          PathCoefficients pc = new PathCoefficients (b, nRays, nOscillators);
          ps.receivers.add(pc);
          return pc.getRxPowerDbm (txPowerDbm);
        }
      }
      PathsSet ps = new PathsSet(a);
      PathCoefficients pc = new PathCoefficients (b, nRays, nOscillators);
      ps.receivers.add(pc);
      paths.add(ps);
      return pc.getRxPowerDbm (txPowerDbm);
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute()
     */
    public double compute() {
      throw new RuntimeException("not supported");
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#compute(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double compute(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      return getRxPower(0, srcRadio, dstRadio);
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.Fading#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location)
     */
    public double computeLimit(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation) {
      throw new RuntimeException("not supported");
    }

  }

} // class: Fading

