001 package org.wdssii.polarmerger;
002
003 import org.apache.commons.logging.Log;
004 import org.apache.commons.logging.LogFactory;
005
006 public class PowerDensityLookup {
007
008 private static Log log = LogFactory.getLog(PowerDensityLookup.class);
009
010 private float beamwidth;
011 private float azimuthalSpacingDegrees;
012 private float gateSpacingKms;
013 private int numGatesSmoothing;
014 private int numBeamsSpread;
015 private float[][][] weights;
016
017 public PowerDensityLookup(float beamwidth, float azimuthalSpacingDegrees, float gatewidthKms, float gateSpacingKms){
018 log.info("Creating PowerDensityLookup");
019 this.beamwidth = beamwidth;
020 this.azimuthalSpacingDegrees = azimuthalSpacingDegrees;
021 this.gateSpacingKms = gateSpacingKms;
022 this.numGatesSmoothing = 1 + Math.round(gatewidthKms / gateSpacingKms);
023 this.numBeamsSpread = 1 + Math.round(beamwidth / azimuthalSpacingDegrees);
024 this.weights = new float[2*numBeamsSpread+1][2*numBeamsSpread+1][2*numGatesSmoothing+1];
025
026 // The lookup stores the weights at various points around a center point
027 float maxPower = 0;
028 float totPower = 0;
029 for (int e=-numBeamsSpread; e <= numBeamsSpread; ++e){
030 float elev_power = computePowerDensity(((float)e)/numBeamsSpread);
031 for (int i=-numBeamsSpread; i <= numBeamsSpread; ++i){
032 float az_power = computePowerDensity(((float)i)/numBeamsSpread);
033 for (int j=-numGatesSmoothing; j <= numGatesSmoothing; ++j){
034 float rd_power = computePowerDensity(((float)j)/numGatesSmoothing);
035 float power = elev_power * az_power * rd_power;
036 weights[e+numBeamsSpread][i+numBeamsSpread][j+numGatesSmoothing] = power;
037 maxPower = Math.max(power, maxPower);
038 totPower += power;
039 }
040 }
041 }
042 for (int i=0; i < weights.length; ++i){
043 for (int j=0; j < weights[0].length; ++j){
044 for (int k=0; k < weights[0][0].length; ++k){
045 weights[i][j][k] /= maxPower;
046 // System.out.println("wt(" + i + "," + j + "," + k + ")=" + weights[i][j][k]);
047 }
048 }
049 }
050
051 log.info("Finished creating PowerDensityLookup");
052 }
053
054 public float getAzimuthalSpacingDegrees() {
055 return azimuthalSpacingDegrees;
056 }
057
058 public float getGateSpacingKms() {
059 return gateSpacingKms;
060 }
061
062
063 public float getBeamwidth() {
064 return beamwidth;
065 }
066
067
068
069 public int getNumBeamsSpread() {
070 return numBeamsSpread;
071 }
072
073
074
075 public int getNumGatesSmoothing() {
076 return numGatesSmoothing;
077 }
078
079
080
081 public float[][][] getWeights() {
082 return weights;
083 }
084
085
086
087 private static float computePowerDensity(float fracOfBeamWidthFromCenter) {
088 // eq. 3.2a of Doviak-Zrnic gives the power density function
089 // theta -- angular distance from the beam axis = frac * beamwidth
090 // sin(theta) is approx theta = dist*beamwidth
091 // lambda (eq.3.2b) = beamwidth*D / 1.27
092 // thus, Dsin(theta)/lambda = D * dist * beamwidth / (D*beamwidth/1.27)
093 // = 1.27 * dist
094 // and the weight is ( J_2(pi*1.27*dist)/(pi*1.27*dist)^2 )^2
095 // dist is between -1/2 and 1/2, so abs(pi*1.27*dist) < 2
096 // in this range, J_2(x) is approximately 1 - exp(-x^2/8.5)
097 double x = (Math.PI * 1.27 * fracOfBeamWidthFromCenter);
098 if (x < 0.0001 && x > -0.0001){
099 x = 0.0001; // avoid divide-by-zero
100 }
101 double wt = (1 - Math.exp(-x * x / 8.5)) / (x * x);
102 wt = wt * wt;
103 return (float) (wt);
104 }
105 }