//////////////////////////////////////////////////
// JIST (Java In Simulation Time) Project
// Timestamp: <PathLoss.java Wed 2004/06/23 09:18:00 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 org.apache.commons.math.MathException;
import org.apache.commons.math.distribution.DistributionFactory;
import org.apache.commons.math.distribution.NormalDistribution;
import org.apache.commons.math.random.AbstractRandomGenerator;
import org.apache.commons.math.random.RandomData;
import org.apache.commons.math.random.RandomDataImpl;
import org.apache.log4j.Logger;

import jist.runtime.JistAPI;
import jist.runtime.Main;
import jist.swans.radio.RadioInfo;
import jist.swans.misc.Location;
import jist.swans.misc.Util;

import jist.swans.Constants;

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

public interface PathLoss
{

  //////////////////////////////////////////////////
  // interface
  //

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

  /**
   * Compute the path loss limit. For deterministic propagation models, it is
   * the same as compute().
   *
   * @param srcRadio source radio information
   * @param srcLocation source location
   * @param dstRadio destination radio information
   * @param dstLocation destination location
   * @return path loss (units: dB)
   */
  double computeLimit(RadioInfo srcRadio, Location srcLocation,
      RadioInfo dstRadio, Location dstLocation);


//  double computeLimit(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo,
//      Location dstLoc, double quantile);


  //////////////////////////////////////////////////
  // interface
  //

  /**
   * Sub class for no path loss.
   */
  public class None implements PathLoss {

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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.PathLoss#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;
    }

//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
//     */
//    public double computeLimit(RadioInfo srcInfo, Location srcLoc,
//        RadioInfo dstInfo, Location dstLoc, double quantile) {
//      return 0;
//    }
  }

  /**
   * Computes free-space path loss. Equivalent to GloMoSim code.
   *
   * @author Rimon Barr &lt;barr+jist@cs.cornell.edu&gt;
   * @since SWANS1.0
   */
  final class FreeSpace implements PathLoss
  {
    // PathLoss interface
    /** {@inheritDoc} */
    public double compute(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation)
    {
      double dist = srcLocation.distance(dstLocation);
      double pathloss = - srcRadio.getGain() - dstRadio.getGain();
      double valueForLog = 4.0 * Math.PI * dist / srcRadio.getWaveLength();
      if (valueForLog > 1.0)
      {
        pathloss += Util.log((float)valueForLog) / Constants.log10 * 20.0;
      }
      return pathloss;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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 compute(srcRadio, srcLocation, dstRadio, dstLocation);
    }

//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
//     */
//    public double computeLimit(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo, Location dstLoc, double quantile) {
//      return compute(srcInfo, srcLoc, dstInfo, dstLoc);
//    }
//
//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#computeCdf(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
//     */
//    public double computeCdf(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo, Location dstLoc, double threshold) {
//      double pathloss = compute(srcInfo, srcLoc, dstInfo, dstLoc);
//
//      return (threshold < -pathloss ? 0. : 1.);
//    }
  } // class: FreeSpace

  /**
   * Computes two-ray path loss. Equivalent to GloMoSim code.
   *
   * @author Rimon Barr &lt;barr+jist@cs.cornell.edu&gt;
   * @since SWANS1.0
   */
  final class TwoRay implements PathLoss
  {
    // PathLoss interface
    /** {@inheritDoc} */
    public double compute(RadioInfo srcRadio, Location srcLocation,
        RadioInfo dstRadio, Location dstLocation)
    {
      double dist = srcLocation.distance(dstLocation);
      double pathloss = - srcRadio.getGain() - dstRadio.getGain();
      double planeEarthLoss = (dist * dist) /
        (srcLocation.getHeight() * dstLocation.getHeight());
      double freeSpaceLoss = 4.0 * Math.PI * dist / srcRadio.getWaveLength();
      if (planeEarthLoss > freeSpaceLoss)
      {
        if (planeEarthLoss > 1.0)
        {
          pathloss += 20.0 * Math.log(planeEarthLoss) / Constants.log10;
        }
      }
      else
      {
        if (freeSpaceLoss > 1.0)
        {
          pathloss += 20.0 * Math.log(freeSpaceLoss) / Constants.log10;
        }
      }
      return pathloss;
    }
    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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 compute(srcRadio, srcLocation, dstRadio, dstLocation);
    }

//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
//     */
//    public double computeLimit(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo, Location dstLoc, double quantile) {
//      return compute(srcInfo, srcLoc, dstInfo, dstLoc);
//    }
//
//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#computeCdf(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
//     */
//    public double computeCdf(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo, Location dstLoc, double threashold) {
//      double pathloss = compute(srcInfo, srcLoc, dstInfo, dstLoc);
//
//      return (threashold < -pathloss ? 0. : 1.);
//    }
  } // class: TwoRay

  // todo: MITRE's pathloss format
  // Time (nearest whole second)    Node A     Node B     Path Loss (dB) Range (meters)
  // End of file is indicated by a -1 in the first column.  (And nothing else on the line.)

  /**
   * Simple log distance propagation modell (without cut-off)
   * This model calculates the reception power with a so-called
   * log-distance propagation model:
   * \f$ L = L_0 + 10 n log_{10}(\frac{d}{d_0})\f$
   *
   * where:
   *  - \f$ n \f$ : the path loss distance exponent
   *  - \f$ d_0 \f$ : reference distance (m)
   *  - \f$ L_0 \f$ : path loss at reference distance (dB)
   *  - \f$ d \f$ : distance (m)
   *  - \f$ L \f$ : path loss (dB)
   *
   * When the path loss is requested at a distance smaller than
   * the reference distance, the tx power is returned.
   *
   * @author kurth
   */
  public static class LogDistance implements PathLoss {

    /** Reference distance in meters (i.e. d0) */
    private double refDist = 1.;

    /** Path loss exponent */
    private double pathlossExp = 2.;

    /**
     * Path loss at the reference point in dB. We treat the value
     * as gain and therefore it's negative. Radio gain is not included
     */
    private double refAttenuationDB;

    private double refWaveLength;

//    private double systemLoss = 1.;

    /**
     * Configure the propagation model. The method recognizes the
     * following DML attributes:
     *
     * @param pathloss_exponent path loss exponent for distance
     *              fading (default: 2 for free space)
     * @param refDist reference distance in meters
     *              (default: 1 meter)
     */
    public LogDistance(double pathloss_exponent,
        double refDist) {
      this.pathlossExp = pathloss_exponent;
      this.refDist = refDist;

      if(this.refDist < 1)
        throw new RuntimeException("DistShadowing: invalid reference_distance "
            + "(must be no less than 1.0)");

      if(this.pathlossExp < 1)
        throw new RuntimeException("DistShadowing: "
          + "invalid pathloss_exponent (must be no less than 1.0)");
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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) {
      double gain = calcPathloss(srcRadio, srcLocation, dstRadio, dstLocation);
      return gain;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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) {
      // Calc path loss
      return calcPathloss(srcRadio, srcLocation, dstRadio, dstLocation);
      }

//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
//     */
//    public double computeLimit(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo,
//        Location dstLoc, double quantile) {
//      return calcPathloss(srcInfo, srcLoc, dstInfo, dstLoc);
//    }

    /**
     * Calculate pathloss a la Friis in dBm.
     *
     * @param srcRadio
     * @param dstRadio
     * @param dist
     * @return pathloss a la Friis in dBm.
     */
    private final double friis(RadioInfo srcRadio, RadioInfo dstRadio, double dist) {
      /*
       * Friis free space equation: (L-system loss)
       *
       *       Pt * Gt * Gr * (lambda^2)
       *   P = --------------------------
       *       (4 * pi * d)^2 * L
       */
      double M = srcRadio.getWaveLength() / (4.0 * Math.PI * dist);
      return ( - 2. * Util.toDB(M) /* systemLoss */);
    }

    /**
     * Calculate the path loss according to the reference distance propagation
     * model.
     *
     * @param srcRadio
     * @param srcLocation
     * @param dstRadio
     * @param dstLocation
     * @return path loss according to the reference distance propagation model
     */
    protected double calcPathloss(RadioInfo srcRadio, Location srcLocation, RadioInfo
        dstRadio, Location dstLocation) {
      // within reference distance, return the transmit power
      double distance = srcLocation.distance(dstLocation);
      if(distance < refDist)
        return 0;

      if (0 == refAttenuationDB) {
        refAttenuationDB = friis(srcRadio, dstRadio, refDist);
        refWaveLength = srcRadio.getWaveLength();
      }

      if (Main.ASSERT)
        Util.assertion(srcRadio.getWaveLength() == refWaveLength);

      return (refAttenuationDB - srcRadio.getGain() - dstRadio.getGain()
          + pathlossExp * Util.toDB(distance/refDist));
    }
  } // class DistShadowing


//  /**
//   * Difference Fading / Shadowing
//   * Log-normal shadowing is a path loss effect. It models the overall signal
//   * attenuation.
//   *
//   * Rician and Rayleigh models refer to the multipath regime. The Rayleigh model
//   * is for diffuse multipath; Rician is for the case where there is a strong
//   * direct ray + diffuse multipath.
//   *
//   * Think of it this way. For the Rayleigh situation, you have a lot of rays of
//   * varying but roughly similar amplitudes adding together as complex phasors.
//   * The real-imaginary probability density function is the 2-D Gaussian with zero
//   * means. The PDF of the envelope of this summation of phasors is a Rayleigh
//   * distribution - basic transformation of random variables.
//   *
//   * Now imagine that one of the rays is much stronger than the others. This will
//   * give the 2-D Gaussian a non-zero mean value, pushing the Gaussian "bump" away
//   * from zero. The PDF of the envelope in this situation is the Rician.
//   *
//   * @author Kurth
//   *
//   * TODO test if computeLimit works!!!!!
//   */
//  public class FriisShadowing implements PathLoss {
//
//    /** logger for field events. */
//    public static final Logger log = Logger.getLogger(FriisShadowing.class.getName());
//
//    /**
//     * path loss exponent
//     */
//    protected double pathlossExp = 2.0;
//
//    /**
//     * shadowing deviation in dB
//     */
//    protected double std_db = 4.0;
//
//    /**
//     * Cache for gaussian.inverseCumulativeProbability(quantile);
//     */
//    protected double quantileProb = .0;
//
//    /**
//     * Used to calculate quantiles.
//     */
//    protected NormalDistribution gaussian;
//
//    /**
//     * Quantile for propagation limit estimation. Must be in [0, 1]. A value
//     * of 0.95 means that 95% of all computations fall within the computed limit.
//     */
//    private double quantile;
//
//
//    public FriisShadowing(double pathlossExp, double std_db, double quantile) {
//      this.pathlossExp = pathlossExp;
//      this.std_db = std_db;
//      this.quantile = quantile;
//
//      DistributionFactory factory = DistributionFactory.newInstance();
//      gaussian = factory.createNormalDistribution();
//
//      try {
//        this.quantileProb = gaussian.inverseCumulativeProbability(quantile);
//      } catch (MathException e) {
//        throw new RuntimeException(e);
//      }
//    }
//
//    /*parent
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#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) {
//
//      // Calc path loss
//      double pathloss = friis(srcRadio, srcLocation, dstRadio, dstLocation);
//
//      // get power loss by adding a log-normal random variable (shadowing)
//      double shadow = Constants.random.nextGaussian() * std_db;
//
//      if (log.isDebugEnabled())
//        log.debug("Src " + srcRadio.getId() + " (" + srcLocation +
//            ") -> Dst " + dstRadio.getId() + " (" + dstLocation +
//            ") has pathloss " + pathloss + " and shadow " + shadow + " = " +
//            pathloss + shadow);
//
//      return (pathloss + shadow);
//    }
//
//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#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) {
//
//      // Calc path loss
//      double pathloss = friis(srcRadio, srcLocation, dstRadio, dstLocation);
//
//      // Estimate the
//      double limit = pathloss - std_db * quantileProb;
//
//      if (log.isDebugEnabled())
//        log.debug("For pathloss " + pathloss + " the " + quantile + " quantile is "
//            + limit);
//
//      return limit;
//      }
//
//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
//     */
//    public double computeLimit(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo, Location dstLoc, double quantile) {
//      // Calc path loss
//      double pathloss = friis(srcInfo, srcLoc, dstInfo, dstLoc);
//
//      // Estimate the
//      double limit = pathloss - std_db * quantileProb;
//
//      if (log.isDebugEnabled())
//        log.debug("For pathloss " + pathloss + " the " + quantile + " quantile is "
//            + limit);
//
//      return limit;
//    }
//
//    /*
//     * (non-Javadoc)
//     * @see jist.swans.field.PathLoss#computeCdf(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
//     */
//    public double computeCdf(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo, Location dstLoc, double pathlossThreashold) {
//      // Calc path loss
//      double pathloss = friis(srcInfo, srcLoc, dstInfo, dstLoc);
//      double cdf = .0;
//
//      try {
//        cdf = gaussian.cumulativeProbability( (pathlossThreashold + pathloss) / std_db);
//      } catch (MathException e) {
//        throw new RuntimeException(e);
//      }
//
//      return cdf;
//    }
//
//    /**
//     * Calculate pathloss a la Friis in dBm.
//     *
//     * @param srcRadio
//     * @param srcLocation
//     * @param dstRadio
//     * @param dstLocation
//     * @return pathloss a la Friis in dBm.
//     */
//    private double friis(RadioInfo srcRadio, Location srcLocation, RadioInfo dstRadio,
//        Location dstLocation) {
//      double dist = srcLocation.distance(dstLocation);
//      /*
//       * Friis free space equation: (L-system loss)
//       *
//       *       Pt * Gt * Gr * (lambda^2)
//       *   P = --------------------------
//       *       (4 * pi * d)^2 * L
//       */
//      double gain = -srcRadio.getGain()
//              - dstRadio.getGain();
//
//      double M = srcRadio.getWaveLength() / (4.0 * Math.PI * dist);
//      return (gain - pathlossExp * Util.toDB(M));
//    }
//
//
//    /*
//     * (non-Javadoc)
//     * @see java.lang.Object#equals(java.lang.Object)
//     */
//    public boolean equals(Object o) {
//      if (this == o) return true;
//      if (o == null || getClass() != o.getClass()) return false;
//
//      final FriisShadowing shadowing = (FriisShadowing) o;
//
//      if (Double.compare(shadowing.pathlossExp, pathlossExp) != 0) return false;
//      if (Double.compare(shadowing.std_db, std_db) != 0) return false;
//
//      return true;
//    }
//
//    /*
//     * (non-Javadoc)
//     * @see java.lang.Object#hashCode()
//     */
//    public int hashCode() {
//      int result;
//      long temp;
//      temp = pathlossExp != +0.0d ? Double.doubleToLongBits(pathlossExp) : 0L;
//      result = (int) (temp ^ (temp >>> 32));
//      temp = std_db != +0.0d ? Double.doubleToLongBits(std_db) : 0L;
//      result = 29 * result + (int) (temp ^ (temp >>> 32));
//      return result;
//    }
//  } // class FriisShadowing
//


  /**
   * Time-Correlated and spatial correlated Shadow Fading
   *
   * This class derives from the Propagation based class to model generic
   * radio signal attenuation due to distance fading as well as the
   * shadowing effect.
   * <pre>Gain_dB = Friis_gain_dB(d0)-10*beta*log10(d/d0)+X_dB</pre>
   * Where Friis_gain_dB(d0) is the signal attenuation (in dB) at
   * reference distance d0 (which is a negative value), beta is the
   * path loss exponent, d>d0 is the distance between the transmitter
   * and the receiver, and X_dB is random variable of Gaussian
   * distribution with zero mean and a standard deviation of delta_dB.
   * The model adds randomness in terms of fluctuation in the radio
   * signal transmission. Typical values of the path loss exponent beta
   * can be set as follows (according to ns-2):
   * <ul>
   *  <li>Outdoor, free space: beta=2
   *  <li>Outdoor, shadowed urban area: beta=2.7~5
   *  <li>Indoor, line of sight: beta=1.6~1.8
   *  <li>Indoor, obstructed: beta=4~6
   * </ul>
   * Typical values of the shadowing deviation delta_dB can be set
   * as follows:
   * <ul>
   *  <li>Outdoor: delta_dB=4~12
   *  <li>Office, hard partition: delta_dB=7
   *  <li>Office, soft partition: delta_dB=9.6
   *  <li>Factory, line of sight: delta_dB=3~6
   *  <li>Factory, obstructed: delta_dB=6.8
   * </ul>
   * The shadowing deviation can be set to be zero in which case we
   * only consider the distance path loss and are free from the
   * probablity part that describes the slow fading.
   *
   * Originated from SWAN (http://www.eg.bucknell.edu/swan/doc/html/index.html)

   * 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).
   *
   *
   * Difference Fading / Shadowing
   * Log-normal shadowing is a path loss effect. It models the overall signal
   * attenuation.
   *
   * Rician and Rayleigh models refer to the multipath regime. The Rayleigh model
   * is for diffuse multipath; Rician is for the case where there is a strong
   * direct ray + diffuse multipath.
   *
   * Think of it this way. For the Rayleigh situation, you have a lot of rays of
   * varying but roughly similar amplitudes adding together as complex phasors.
   * The real-imaginary probability density function is the 2-D Gaussian with zero
   * means. The PDF of the envelope of this summation of phasors is a Rayleigh
   * distribution - basic transformation of random variables.
   *
   * Now imagine that one of the rays is much stronger than the others. This will
   * give the 2-D Gaussian a non-zero mean value, pushing the Gaussian "bump" away
   * from zero. The PDF of the envelope in this situation is the Rician.
   *
   * @author Kurth
   *
   * TODO spatial correlation
   */
  public static class Shadowing implements PathLoss {

    /** logger for field events. */
    private static final Logger log = Logger.getLogger(Shadowing.class.getName());

    /**
     *
     * Used by Rayleigh and Ricean fading and shadow fading: two nodes
     * should be correlated for a short time period.
     */
    public static class CorrelatedFading {
      public CorrelatedFading(double fading, long expires) {
        this.fading = fading;
        this.expires = expires;
      }
      double fading;
      long expires;
    }

    /** Shadowing deviation */
    protected double shadowingStdev = .0;

    /** Beyond this time, we assume slow fading is independent. */
    protected long shadowingCoherenceTime = 1000 * Constants.MILLI_SECOND;

    /** whether shadow intervals should be exponential distributed */
    protected boolean exponential = true;

    /** Whether links and link correlation should be symmetric */
    protected boolean symmetric = true;

    /**
     * Quantile for propagation limit estimation. Must be in [0, 1]. A value
     * of 0.95 means that 95% of all computations fall within the computed limit.
     */
    protected double quantile;

    /** deterministic pathloss like free-space or two-ray */
    protected PathLoss pathloss;


    // internal

    /** Storing half of the matrix in an array: indexed by node ids. */
    private CorrelatedFading[][] shadowingCache = null;

    /**
     * Cache for gaussian.inverseCumulativeProbability(quantile);
     */
    private double quantileProb = .0;

    /** distribution of shadowing intervals */
    private RandomData shadowIntervalDist;

    /**
     * Configure the propagation model. The method recognizes the
     * following DML attributes:
     *
     * @param shadowingStdev standard deviation for the Gaussian
     *              shadow fading in dB (default: 0 for no shadowing)
     * @param shadowingCoherenceTime duration within which slow
     *              fading is treated as correlated (in milliseconds, default=1000).
     * @param quantile
     */
    public Shadowing(PathLoss pathloss, double shadowingStdev,
        long shadowingCoherenceTime, boolean exponential, boolean symmetric,
        double quantile, int noNodes) {
      this.pathloss = pathloss;

      this.shadowingStdev = shadowingStdev;
      this.shadowingCoherenceTime = shadowingCoherenceTime;

      if(this.shadowingStdev < 0)
        throw new RuntimeException("DistShadowing: "
            + "invalid shadowing_stdev (must be positive)");

      if(this.shadowingCoherenceTime < 0)
        throw new RuntimeException("DistShadowing: "
            + "invalid slow_fading_coherence_time (must be positive)");

      resizeCache(noNodes);
      this.exponential = exponential;
      this.symmetric = symmetric;
      this.quantile = quantile;

      if (this.exponential)
        shadowIntervalDist = new RandomDataImpl(new AbstractRandomGenerator() {
          public double nextDouble() {
            return Constants.random.nextDouble();
          }
          public void setSeed(long seed) {
            Constants.random.setSeed(seed);
          }
        });

      DistributionFactory factory = DistributionFactory.newInstance();
      NormalDistribution gaussian = factory.createNormalDistribution();

      try {
        this.quantileProb = gaussian.inverseCumulativeProbability(quantile);
      } catch (MathException e) {
        throw new RuntimeException(e);
      }
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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) {
      double gain = pathloss.compute(srcRadio, srcLocation, dstRadio, dstLocation);
      gain += getShadowDB(srcRadio, dstRadio);
      if(gain < 0)
        gain = 0; // make sure it's not negative because of randomness.

      //gain += fast_fading_db(fromidx, toidx);
      return gain;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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 computeLimit(srcRadio, srcLocation, dstRadio, dstLocation, quantile);
    }


    public double computeLimit(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo, Location dstLoc, double quantile) {
      // Calc path loss
      double pathloss = this.pathloss.compute(srcInfo, srcLoc, dstInfo, dstLoc);

      // Estimate the
      double limit = pathloss - shadowingStdev * quantileProb;

      if (log.isDebugEnabled())
        log.debug("For pathloss " + pathloss + " the " + quantile + " quantile is "
            + limit);

      return limit;
    }

    /**
     * Resize the shadowing cache.
     *
     * @param idxMax the new size of the cache.
     */
    private void resizeCache(int idxMax) {
      if (Main.ASSERT)
        Util.assertion(null == shadowingCache || idxMax >= shadowingCache.length);

      CorrelatedFading[][] cache = new CorrelatedFading[idxMax+1][idxMax+1];
      for (int x=0; null != shadowingCache && x < shadowingCache.length; x++) {
        for (int y=0; y < shadowingCache[x].length; y++) {
          cache[x][y] = shadowingCache[x][y];
        }
      }
      shadowingCache = cache;
    }

    /**
     * Compute correlated shadowing. The link from src to dst is assumed to be
     * correlated in time for the interval slow_fading_cache.
     *
     * @param srcRadio
     * @param dstRadio
     * @return correlated shadowing in dB.
     */
    private double getShadowDB(RadioInfo srcRadio, RadioInfo dstRadio) {
      if(shadowingStdev == 0)
        return 0;

      int n1 = srcRadio.getId();
      int n2 = dstRadio.getId();
      if(this.symmetric && n1 > n2) {
        int tmp = n1; n1 = n2; n2 = tmp;
      }

      if(n1 >= shadowingCache.length || n2 >= shadowingCache[n1].length)
        resizeCache(Math.max(n1, n2));

      if(null == shadowingCache[n1][n2]) {
        shadowingCache[n1][n2] = new CorrelatedFading(shadowingStdev *
            Constants.random.nextGaussian(), JistAPI.getTime() +
            (!exponential ? shadowingCoherenceTime :
                (long)shadowIntervalDist.nextExponential(shadowingCoherenceTime)));

        // Prevents from re-calculating
        if (shadowingCoherenceTime == JistAPI.END)
          shadowingCache[n1][n2].expires = JistAPI.END;
      }
      else if(.0 == shadowingCoherenceTime
          ||JistAPI.getTime() > shadowingCache[n1][n2].expires) {
        shadowingCache[n1][n2].fading = shadowingStdev * Constants.random.nextGaussian();
        shadowingCache[n1][n2].expires = JistAPI.getTime() +
          (!exponential ? shadowingCoherenceTime :
            (long)shadowIntervalDist.nextExponential(shadowingCoherenceTime));
      }

      return shadowingCache[n1][n2].fading;
    }

  } // class DistShadowing


  /**
   * @deprecated use {@link PathLoss.Shadowing}
   */
  public static class DistShadowing implements PathLoss {

    /** logger for field events. */
    public static final Logger log = Logger.getLogger(DistShadowing.class.getName());

    /**
     *
     * Used by Rayleigh and Ricean fading and shadow fading: two nodes
     * should be correlated for a short time period.
     */
    public static class CorrelatedFading {
      public CorrelatedFading(double fading, long recorded_time) {
        this.fading = fading;
        this.recorded_time = recorded_time;
      }
      double fading;
      long recorded_time;
    }

    /** Reference distance in meters (i.e. d0) */
    private double refDist = 1.;

    /** Path loss exponent */
    private double pathlossExp = 2.;

    /** Shadowing deviation */
    private double shadowingStdev = .0;

//    private double systemLoss = 1.;

    /**
     * Path loss at the reference point in dB. We treat the value
     * as gain and therefore it's negative.
     */
    private double refAttenuationDB;

    /** Beyond this time, we assume slow fading is independent. */
    private long shadowingCoherenceTime = 1000 * Constants.MILLI_SECOND;

    /** Whether links and link correlation should be symmetric */
    private boolean symmetric = true;

    /** Storing half of the matrix in an array: indexed by node ids. */
    private CorrelatedFading[][] shadowingCache = null;

    /**
     * Cache for gaussian.inverseCumulativeProbability(quantile);
     */
    protected double quantileProb = .0;

    /**
     * Used to calculate quantiles.
     */
    protected NormalDistribution gaussian;

    /**
     * Quantile for propagation limit estimation. Must be in [0, 1]. A value
     * of 0.95 means that 95% of all computations fall within the computed limit.
     */
    private double quantile;

    /**
     * Configure the propagation model. The method recognizes the
     * following DML attributes:
     *
     * @param pathloss_exponent path loss exponent for distance
     *              fading (default: 2 for free space)
     * @param shadowing_stdev standard deviation for the Gaussian
     *              shadow fading in dB (default: 0 for no shadowing)
     * @param reference_distance reference distance in meters
     *              (default: 1 meter)
     * @param slow_fading_coherence_time duration within which slow
     *              fading is treated as correlated (in milliseconds, default=1000).
     * @param quantile
     */
    public DistShadowing(double pathloss_exponent, double shadowing_stdev,
        double reference_distance, long slow_fading_coherence_time,
        int noNodes, boolean symmetric, double quantile) {
      this.pathlossExp = pathloss_exponent;
      this.shadowingStdev = shadowing_stdev;
      this.refDist = reference_distance;
      this.shadowingCoherenceTime = slow_fading_coherence_time;

      if(this.refDist < 1)
        throw new RuntimeException("DistShadowing: invalid reference_distance "
            + "(must be no less than 1.0)");

      if(this.pathlossExp < 1)
        throw new RuntimeException("DistShadowing: "
          + "invalid pathloss_exponent (must be no less than 1.0)");

      if(this.shadowingStdev < 0)
        throw new RuntimeException("DistShadowing: "
            + "invalid shadowing_stdev (must be positive)");

      if(this.shadowingCoherenceTime < 0)
        throw new RuntimeException("DistShadowing: "
            + "invalid slow_fading_coherence_time (must be positive)");

      resizeCache(noNodes);
      this.symmetric = symmetric;
      this.quantile = quantile;

      DistributionFactory factory = DistributionFactory.newInstance();
      gaussian = factory.createNormalDistribution();

      try {
        this.quantileProb = gaussian.inverseCumulativeProbability(quantile);
      } catch (MathException e) {
        throw new RuntimeException(e);
      }
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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) {
      double gain = calcPathloss(srcRadio, srcLocation, dstRadio, dstLocation);
      gain += getShadowDB(srcRadio, dstRadio);
      if(gain < 0)
        gain = 0; // make sure it's not negative because of randomness.

      //gain += fast_fading_db(fromidx, toidx);
      return gain;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#computeCdf(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
     */
    public double computeCdf(RadioInfo srcInfo, Location srcLoc, RadioInfo
        dstInfo, Location dstLoc, double pathlossThreashold) {
      // Calc path loss
      double pathloss = calcPathloss(srcInfo, srcLoc, dstInfo, dstLoc);
      double cdf = .0;

      try {
        cdf = gaussian.cumulativeProbability( (pathlossThreashold + pathloss) /
            shadowingStdev);
      } catch (MathException e) {
        throw new RuntimeException(e);
      }

      return cdf;
    }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#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) {

      // Calc path loss
      double pathloss = calcPathloss(srcRadio, srcLocation, dstRadio, dstLocation);

      // Estimate the
      double limit = pathloss - shadowingStdev * quantileProb;

      if (log.isDebugEnabled())
        log.debug("For pathloss " + pathloss + " the " + quantile + " quantile is "
            + limit);

      return limit;
      }

    /*
     * (non-Javadoc)
     * @see jist.swans.field.PathLoss#computeLimit(jist.swans.radio.RadioInfo, jist.swans.misc.Location, jist.swans.radio.RadioInfo, jist.swans.misc.Location, double)
     */
    public double computeLimit(RadioInfo srcInfo, Location srcLoc, RadioInfo dstInfo, Location dstLoc, double quantile) {
      // Calc path loss
      double pathloss = calcPathloss(srcInfo, srcLoc, dstInfo, dstLoc);

      // Estimate the
      double limit = pathloss - shadowingStdev * quantileProb;

      if (log.isDebugEnabled())
        log.debug("For pathloss " + pathloss + " the " + quantile + " quantile is "
            + limit);

      return limit;
    }

    /**
     * Resize the shadowing cache.
     *
     * @param idxMax the new size of the cache.
     */
    private void resizeCache(int idxMax) {
      if (Main.ASSERT)
        Util.assertion(null == shadowingCache || idxMax >= shadowingCache.length);

      CorrelatedFading[][] cache = new CorrelatedFading[idxMax+1][idxMax+1];
      for (int x=0; null != shadowingCache && x < shadowingCache.length; x++) {
        for (int y=0; y < shadowingCache[x].length; y++) {
          cache[x][y] = shadowingCache[x][y];
        }
      }
      shadowingCache = cache;
    }

    /**
     * Compute correlated shadowing. The link from src to dst is assumed to be
     * correlated in time for the interval slow_fading_cache.
     *
     * @param srcRadio
     * @param dstRadio
     * @return correlated shadowing in dB.
     */
    private double getShadowDB(RadioInfo srcRadio, RadioInfo dstRadio) {
      if(shadowingStdev == 0)
        return 0;

      int n1 = srcRadio.getId();
      int n2 = dstRadio.getId();
      if(this.symmetric && n1 > n2) {
        int tmp = n1; n1 = n2; n2 = tmp;
      }

      if(n1 >= shadowingCache.length || n2 >= shadowingCache[n1].length)
        resizeCache(Math.max(n1, n2));

      if(null == shadowingCache[n1][n2]) {
        shadowingCache[n1][n2] = new CorrelatedFading(shadowingStdev *
            Constants.random.nextGaussian(), JistAPI.getTime());
      }
      else if(.0 == shadowingCoherenceTime
          ||JistAPI.getTime() > shadowingCache[n1][n2].recorded_time + shadowingCoherenceTime) {
        shadowingCache[n1][n2].fading = shadowingStdev * Constants.random.nextGaussian();
        shadowingCache[n1][n2].recorded_time = JistAPI.getTime();
      }

      return shadowingCache[n1][n2].fading;
    }

    /**
     * Calculate pathloss a la Friis in dBm.
     *
     * @param srcRadio
     * @param dstRadio
     * @param dist
     * @return pathloss a la Friis in dBm.
     */
    private double friis(RadioInfo srcRadio, RadioInfo dstRadio, double dist) {
      /*
       * Friis free space equation: (L-system loss)
       *
       *       Pt * Gt * Gr * (lambda^2)
       *   P = --------------------------
       *       (4 * pi * d)^2 * L
       */
      double gain = -srcRadio.getGain() - dstRadio.getGain();

      double M = srcRadio.getWaveLength() / (4.0 * Math.PI * dist);
      return (gain - 2. * Util.toDB(M) /* systemLoss */);
    }

    /**
     * Calculate the path loss according to the reference distance propagation
     * model.
     *
     * @param srcRadio
     * @param srcLocation
     * @param dstRadio
     * @param dstLocation
     * @return path loss according to the reference distance propagation model
     */
    private double calcPathloss(RadioInfo srcRadio, Location srcLocation, RadioInfo
        dstRadio, Location dstLocation) {
      // Calc path loss
      refAttenuationDB = friis(srcRadio, dstRadio, refDist);

      double distance = srcLocation.distance(dstLocation);
      if(distance < 0.5)
        distance = 0.5;

      return (refAttenuationDB + pathlossExp * Util.toDB(distance/refDist));
    }

  } // class DistShadowing

} // class: PathLoss

