ttlayer.java

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

                           (thickness[idx]*sineOfAngle)/Math.sqrt(vSquared[srcLayerIdx]/vSquared[idx] - sineOfAngle*sineOfAngle);
                }

                double rangeErr  = range - estRange;
                if (Math.abs(rangeErr) <= 0.02) {
                    return sineOfAngle;
                }

                if (rangeErr < 0.) {
                    XBIG = XTR;
                    deltaBig = estRange;
                }
                else {
                    XLIT = XTR;
                    deltaSmall = estRange;
                }
            } // end of iteration loop

            XTR = 0.5 * (XBIG + XLIT);
            sineOfAngle = XTR/Math.sqrt(XTR*XTR + srcLayerZSquared);
            return sineOfAngle;
    }

//  Travel time & derivatives for direct wave below first layer
    private void calcDirectWaveResults(double sineOfAngle) {
        double directTT = srcLayerZ/(velocities[srcLayerIdx] * Math.sqrt(1.0 - sineOfAngle*sineOfAngle));
        for (int idx = 0; idx < srcLayerIdx; idx ++) {
            directTT += (thickness[idx] * velocities[srcLayerIdx]) /
                        (vSquared[idx] * Math.sqrt(vSquared[srcLayerIdx]/vSquared[idx] - sineOfAngle*sineOfAngle) );
        }
        if ( (srcLayerIdx == maxLayerIdx) || (directTT < minHeadWaveTravelTime) ) {
            double SRR = Math.sqrt(1. - sineOfAngle*sineOfAngle);
            double SRT = SRR * SRR * SRR;
            double ALPHA = srcLayerZ/SRT;
            double BETA = srcLayerZ * sineOfAngle/(velocities[srcLayerIdx] * SRT);

            for (int idx = 0; idx < srcLayerIdx; idx++) {
                double STK = Math.sqrt(vSquared[srcLayerIdx]/vSquared[idx] - sineOfAngle*sineOfAngle);
                STK = STK * STK * STK;
                double VTK = thickness[idx] / ( vSquared[idx] * STK );
                ALPHA += VTK * vSquared[srcLayerIdx];
                BETA += VTK * velocities[srcLayerIdx] * sineOfAngle;
            }

            calcTT = directTT;
            dTdR = BETA/ALPHA;
            dTdZ = (1.0 - velocities[srcLayerIdx] * sineOfAngle * dTdR)/(velocities[srcLayerIdx] * SRR);
            incAngleDeg = toAngleOfIncidenceDegrees(sineOfAngle);
        }
        else {
            setHeadWaveResults(fastestLayerIdx);
        }
        return;
    }


    private void setHeadWaveResults(int layerIdx) {
        calcTT = layerHeadWaveTravelTime[layerIdx];
        dTdR = 1.0/velocities[layerIdx];
        dTdZ =  - Math.sqrt(vSquared[layerIdx] - vSquared[srcLayerIdx])/(velocities[layerIdx] * velocities[srcLayerIdx]);
        incAngleDeg = toAngleOfIncidenceDegrees(- velocities[srcLayerIdx]/velocities[layerIdx]);
    }

    private double toAngleOfIncidenceDegrees(double sineOfAngle) {
        if ( Math.abs(sineOfAngle) > 1.0 ) sineOfAngle = 1.0;
        double angleDeg = Math.toDegrees(Math.asin(sineOfAngle)); // * 57.29578
        if (angleDeg < 0.) angleDeg += 180.;
        angleDeg = 180. - angleDeg;
        return angleDeg;
    }

/**
*  Setup travel time matrices, matrices are layerUpDownPathTime[srcIdx][headIdx] AND layerUpDownHorizDist[srcIdx][headIdx].
*  Source is at top of srcIdx, ray is critically refracted at top of headIdx (with velocity=velocities[headIdx]
*  layerUpDownPathTime is travel time for that path, layerUpDownHorizDist is surface distance traveled from source along path.
*/
    private void initModelCalculationBuffers() {

        createModelCalculationBuffers();

        // Cache velocities squared for calculations
        for (int idx = 0; idx < nLayers; idx++) {
            vSquared[idx] = velocities[idx] * velocities[idx];
        }

        // Save layer thicknesses for calculations
        int maxIdx = nLayers - 1;
        for (int idx = 0; idx < maxIdx; idx++) {
            thickness[idx] = depths[idx+1] - depths[idx];
        }

        initModelHeadWaveUpDownPathDataBuffers() ;
    }

    // Compute layerUpDownPathTime and layerUpDownHorizDist for src at top of source layer=srcIdx
    private void initModelHeadWaveUpDownPathDataBuffers() {
        // For each src idx set the up-down path leg contribution 
        for (int srcIdx = 0; srcIdx < nLayers; srcIdx++) {
            // Ray refracted along underlying layer identified by headIdx
            for (int headIdx = 0; headIdx < nLayers; headIdx++) {
                pathLegsThruLayer[srcIdx][headIdx] = 1.0;                        // headIdx <  srcIdx up leg only
                if (headIdx >= srcIdx) pathLegsThruLayer[srcIdx][headIdx] = 2.0; // headIdx >= srcIdx up & down legs
            }
        }

        // Zero-out traveltime and distance calculation buffers for current velocity model.
        for (int srcIdx = 0; srcIdx < nLayers; srcIdx++) {
            for (int headIdx = 0; headIdx < nLayers; headIdx++) {
                layerUpDownPathTime[srcIdx][headIdx] = 0.;
                layerUpDownHorizDist[srcIdx][headIdx] = 0.;
            }
        }

        // Src layer reference index
        for (int srcIdx = 0; srcIdx < nLayers; srcIdx++) {
            HeadWaveLayerLoop:
            for (int headIdx = srcIdx; headIdx < nLayers; headIdx++) {  // do underlying layer headwaves
                if (headIdx == 0) continue HeadWaveLayerLoop;
                int lastIdx = headIdx; 
                for (int idx = 0; idx < lastIdx; idx++) {  // calc tt and dist from src thru layers above refracting layer.
                    double velocityDiffSqrt = Math.sqrt(vSquared[headIdx] - vSquared[idx]);
                    double time = thickness[idx] * velocityDiffSqrt/(velocities[idx] * velocities[headIdx]);
                    double dist = thickness[idx] * velocities[idx]/velocityDiffSqrt;
                    layerUpDownPathTime[srcIdx][headIdx] += pathLegsThruLayer[srcIdx][idx] * time;
                    layerUpDownHorizDist[srcIdx][headIdx] += pathLegsThruLayer[srcIdx][idx] * dist;
                } // end of lastIdx loop
            } // end of headIdx inner loop
        } // end of outer srcIdx loop

    } // end of method

    /** Return travel-time (secs) through this model for a receiver at distance
     * 'range' in km. Note that adjustment for receiver elevation is not possible. */
    public double pTravelTime(double range) {
        if (debug) System.out.println("\n****************pTravelTime*********************");
        calculateResults(range);
        return calcTT;
    }

    public double getDtDz() { return dTdZ; }
    public double getDtDr() { return dTdR; }

    public double getSlowness() { return dTdR; }
    public double getAngleOfIncidence() { return incAngleDeg; }

    public static void outputSoCalTable( double startRangeKm,    double rangeInc, int rangeElements,
                                    double startSrcDepthKm, double depthInc, int depthElements) {
        outputTable( UniformFlatLayerVelocityModel.SOCAL_DEFAULT,
                     startRangeKm,    rangeInc, rangeElements,
                     startSrcDepthKm, depthInc, depthElements) ;
    }

    public static void outputTable(UniformFlatLayerVelocityModelIF vm,
                               double startRangeKm,    double rangeInc, int rangeElements,
                               double startSrcDepthKm, double depthInc, int depthElements) {
        TTLayer ttlayer = new TTLayer(vm);
        ttlayer.outputTable(startRangeKm, rangeInc, rangeElements,
                            startSrcDepthKm, depthInc, depthElements);
    }

    protected void outputTable(double startRangeKm, double rangeInc, int rangeElements,
              double startSrcDepthKm, double depthInc, int depthElements) {
        double [] depthKm = new double[depthElements];
        double depth = startSrcDepthKm;
        for (int depthIdx = 0; depthIdx < depthElements; depthIdx++) {
            depthKm[depthIdx] = depth;
            depth += depthInc;
        }
        StringBuffer depthStringBuffer = new StringBuffer(132);
        depthStringBuffer.append("  V  ");
        Concat cc  = new Concat();
        for (int depthIdx = 0; depthIdx < depthElements; depthIdx++) {
            cc.format(depthStringBuffer,depthKm[depthIdx], 3, 2);
//          depthStringBuffer.append(" ");
        }
        double dist = startRangeKm;
        double [] range = new double [rangeElements];
        for (int rangeIdx = 0; rangeIdx < rangeElements; rangeIdx++) {
             range[rangeIdx] = dist;
             dist += rangeInc;
        }
        
        double [][] tt = new double [depthElements][rangeElements];
        for (int depthIdx = 0; depthIdx < depthElements; depthIdx++) {
            this.setSource(new LatLonZ(0.,0., depthKm[depthIdx]));
            for (int rangeIdx = 0; rangeIdx < rangeElements; rangeIdx++) {
                tt[depthIdx][rangeIdx] = this.pTravelTime(range[rangeIdx]);
            }
        }
        if (! debug) {
            System.out.println("Range");
            System.out.println("  |  Depth ->");
            System.out.println(depthStringBuffer.toString());
            for (int rangeIdx = 0; rangeIdx < rangeElements; rangeIdx++) {
                StringBuffer sb = new StringBuffer(132);
//                cc.format(sb, (long) range[rangeIdx],3);
                cc.format(sb, range[rangeIdx], 3, 1);
                for (int depthIdx = 0; depthIdx < depthElements; depthIdx++) {
                    cc.format(sb,tt[depthIdx][rangeIdx],3,2);
//                  sb.append(" ");
                }
                System.out.println(sb.toString());
            }
        }
    }

    public static void main(String [] args) {
        TTLayer ttlayer = new TTLayer(UniformFlatLayerVelocityModel.SOCAL_DEFAULT);
        double depthKm = 0.;
        int depthElements = args.length;
        if (depthElements > 0) {
           depthKm = Double.parseDouble(args[0]);
            ttlayer.outputTable(0., 5., 51, depthKm, 0., 1);
        }
        else {
            ttlayer.outputTable(0., 5., 51, 0., 2., 18);
        }
/*
        LatLonZ receiver = new LatLonZ(34.2, -121., 0.);
        range = GeoidalConvert.horizontalDistanceKmBetween(ttlayer.getSource(), receiver);
        tt = ttlayer.pTravelTime(receiver);
        System.out.println("LatLonZ Range: " + range + " time: " + tt);
*/
    }

} // end of class