ttlayer.java

一个用java写的地震分析软件(无源码)

package org.trinet.util.velocitymodel;
import org.trinet.util.gazetteer.*;

/**
* Derived from TTLAYR J.P. Eaton's FORTRAN algoritim for computing traveltime and derivatives
* for a uniform flat layer crustal model.
* The source depth is in km with a positive value down. Receiver depth is ALWAYS
* considered to be 0.0km.
* Class main(String []) parses the input line a source depth
* for which to calculate traveltime to ranges of 0 to 250 km.<p>
* The main method uses the HadleyKanamori Southern California flat layer velocity model.
* The ranges and calculated travel times are output to System.out.
*/
public class TTLayer extends UniformFlatLayerTravelTimeGenerator {
    private static boolean debug = false;
    private static final int FIRST_LAYER_IDX = 0;

// value derived from input
    private double range;
    private double minHeadWaveTravelTime = 99999.99;
    private double maxCriticalAngleRange;

// derived values calculated for return
    private double calcTT;
    private double dTdZ;
    private double dTdR;
    private double incAngleDeg;

// buffer loop index bounds
    private int nLayers;
    private int maxLayerIdx;
    private int srcLayerIdx;
    private int layerBelowSrcIdx;
    private int fastestLayerIdx;

// source depths (down is positive)
    private double srcZ;
    private double srcLayerZ;
    private double srcLayerZSquared;

// buffers for calculations
    private double []   srcZCorrLayerUpDownPathTime;
    private double []   srcZCorrLayerUpDownHorizDist;

    private double []   layerHeadWaveTravelTime;
    private double []   vSquared;
    private double []   thickness;
    private double [][] pathLegsThruLayer;
    private double [][] layerUpDownPathTime;
    private double [][] layerUpDownHorizDist;

// set by constructor input parameters
    private double [] velocities;
    private double [] depths;

    public TTLayer(UniformFlatLayerVelocityModelIF vm) {
        initModel(vm);  // must initModel first
    }

    public TTLayer(UniformFlatLayerVelocityModelIF vm, Geoidal source) {
        initModel(vm);  // order dependent, must initModel first
        setSource(source);
    }

    private void initModel(UniformFlatLayerVelocityModelIF vm) {
        super.setModel(vm);
        velocities = vm.getVelocities();
        depths = vm.getBoundaryDepths();
        initModelCalculationBuffers();
    }

    public void setModel(VelocityModelIF vm) {
        initModel((UniformFlatLayerVelocityModelIF) vm);
        if (getSource() != null) initCalcDataForSource();
    }

    private void createModelCalculationBuffers() {
        this.nLayers = ((UniformFlatLayerVelocityModelIF) getModel()).getLayerCount();
        this.maxLayerIdx = nLayers - 1;
        srcZCorrLayerUpDownPathTime = new double [nLayers];
        srcZCorrLayerUpDownHorizDist = new double [nLayers];
        layerHeadWaveTravelTime = new double [nLayers];
        vSquared = new double [nLayers];
        layerUpDownPathTime = new double [nLayers][nLayers];
        layerUpDownHorizDist  = new double [nLayers][nLayers];
        thickness = new double [nLayers];
        pathLegsThruLayer = new double [nLayers][nLayers]; 
    }

    public void setSource(Geoidal source) {
        super.setSource(source);
        if (getModel() == null) throw new NullPointerException("VelocityModel null, invoke setVelocity() first.");
        initCalcDataForSource();
    }

    private void initCalcDataForSource() {
        srcZ = Math.abs(getSource().getZ());    // positive depth downwards into model
        srcLayerIdx = getSourceLayerIndex();

        srcLayerZ = srcZ - depths[srcLayerIdx]; // depth of source in source layer: srcLayerZ 
        srcLayerZSquared = srcLayerZ*srcLayerZ + 0.000001;

        layerBelowSrcIdx = srcLayerIdx;
        if (srcLayerIdx <  maxLayerIdx) {      // srcZ < maxLayerZ: get maxCriticalAngleRange, minRefractTime, minRefractDist
            layerBelowSrcIdx++;
            if (getModel() == null) throw new NullPointerException("VelocityModel must be set not null before setting source.");
            initSourceDepthCorrectedUpDownPathData();
        }
    }

    private int getSourceLayerIndex() {
        int index = nLayers - 1;
        for (int idx = 0; idx < nLayers; idx++) {
            if (depths[idx] > srcZ) {
                index = idx-1;
                break;
            }
        }
        return index;
    }

    private void initSourceDepthCorrectedUpDownPathData() {
        for (int idx = layerBelowSrcIdx; idx < nLayers; idx++) {
                double velocityDiffSqrt = Math.sqrt(vSquared[idx] - vSquared[srcLayerIdx]);
                srcZCorrLayerUpDownPathTime[idx] =
                    layerUpDownPathTime[srcLayerIdx][idx] - srcLayerZ*velocityDiffSqrt/(velocities[idx]*velocities[srcLayerIdx]);
                srcZCorrLayerUpDownHorizDist[idx] =
                    layerUpDownHorizDist[srcLayerIdx][idx] - srcLayerZ*velocities[srcLayerIdx]/velocityDiffSqrt;
        }
        maxCriticalAngleRange = ( velocities[layerBelowSrcIdx] * velocities[srcLayerIdx] ) *
                 ( srcZCorrLayerUpDownPathTime[layerBelowSrcIdx] - layerUpDownPathTime[srcLayerIdx][srcLayerIdx] ) /
                 ( velocities[layerBelowSrcIdx] - velocities[srcLayerIdx] ) ;
    }

//  Direct wave if focus in last layer find smallest possible refracted arrival time 
    private void calculateResults(double distanceFromSource) {
        this.range = distanceFromSource;
        fastestLayerIdx = findFastestHeadWaveLayerIndex(range);
        if ( (srcLayerIdx == maxLayerIdx) || (range < maxCriticalAngleRange)) {
            if (srcLayerIdx == FIRST_LAYER_IDX) doFirstLayerSrc();
            else doDirectWave();
        }
        else {
            setHeadWaveResults(fastestLayerIdx);
        }
        return;
    }

    private int findFastestHeadWaveLayerIndex(double range) {
        for (int idx = layerBelowSrcIdx; idx < nLayers; idx++) {
            layerHeadWaveTravelTime[idx]  =  srcZCorrLayerUpDownPathTime[idx] + range/velocities[idx];
        }

        int fastestLayerIdx = 0;
        minHeadWaveTravelTime = 99999.99;
        for ( int idx = layerBelowSrcIdx; idx < nLayers; idx++) {
            if ( (layerHeadWaveTravelTime[idx] > minHeadWaveTravelTime) || (srcZCorrLayerUpDownHorizDist[idx] > range) ) continue;
            fastestLayerIdx = idx;
            minHeadWaveTravelTime = layerHeadWaveTravelTime[idx];
        }
        return fastestLayerIdx;
    }

//  Calculation for a direct wave from source located in the first layer 
    private void doFirstLayerSrc() {
        double dist = Math.sqrt(srcZ*srcZ + range*range);
        double directTT = dist/velocities[FIRST_LAYER_IDX];
        if (directTT >= minHeadWaveTravelTime) {
            setHeadWaveResults(fastestLayerIdx);   // set results here
        }
        calcTT = directTT; // set result here 
        dTdR = range/(velocities[FIRST_LAYER_IDX] * dist);
        dTdZ = srcZ/(velocities[FIRST_LAYER_IDX] * dist);
        incAngleDeg = toAngleOfIncidenceDegrees(range/dist);
    }

//  Find a direct wave that WILL emerge at the current range
    private void doDirectWave() {
        double XBIG = range;
        double XLIT = range * srcLayerZ/srcZ;
        double UB = XBIG/Math.sqrt(XBIG*XBIG + srcLayerZSquared);
        double UL = XLIT/Math.sqrt(XLIT*XLIT + srcLayerZSquared);

        double deltaBig = srcLayerZ * UB/Math.sqrt(1.000001 - UB*UB);
        double deltaSmall = srcLayerZ * UL/Math.sqrt(1.000001 - UL*UL);

        for (int idx = 0; idx < srcLayerIdx; idx++) {
            deltaBig   +=  (thickness[idx] * UB)/Math.sqrt(vSquared[srcLayerIdx]/vSquared[idx] - UB*UB);
            deltaSmall +=  (thickness[idx] * UL)/Math.sqrt(vSquared[srcLayerIdx]/vSquared[idx] - UL*UL);
        }

        double sineOfAngle = extrapolateSineOfAngleForDirectWave(XBIG, XLIT, deltaBig, deltaSmall);
        if (1.0 - sineOfAngle > 0.0002) {
            calcDirectWaveResults(sineOfAngle);
            return;
        }

        //  If sineOfAngle is close to 1.0, compute traveltime for head wave along the top of layer srcLayerIdx
        double tt = layerUpDownPathTime[srcLayerIdx][srcLayerIdx] + range/velocities[srcLayerIdx];

        if ( (srcLayerIdx == maxLayerIdx) || (tt < minHeadWaveTravelTime) ) { // results for head wave at top of src layer
            calcTT = tt;
            dTdR = 1.0/velocities[srcLayerIdx];
            dTdZ = 0.;
            incAngleDeg = toAngleOfIncidenceDegrees(0.9999999);
        }
        else {
            setHeadWaveResults(fastestLayerIdx);  // use fastest head wave traveltime
        }
        return; 
    }

    private double extrapolateSineOfAngleForDirectWave(double XBIG, double XLIT, double deltaBig, double deltaSmall) {
            double XTR = 0.;
            double sineOfAngle = 0.;
            for (int iteration = 1; iteration <= 26; iteration++) {

                if ( (iteration > 10) && ( (1.0 - sineOfAngle) < 0.0002) ) {
                     return sineOfAngle;
                }
                if (deltaBig - deltaSmall < 0.02) break;

                XTR = XLIT + (range - deltaSmall) * (XBIG - XLIT)/(deltaBig - deltaSmall);
                sineOfAngle = XTR/Math.sqrt(XTR*XTR + srcLayerZSquared);

                double estRange = srcLayerZ * sineOfAngle/Math.sqrt(1.000001 - sineOfAngle*sineOfAngle);
                for (int idx = 0; idx < srcLayerIdx; idx++) {
                    estRange +=