| 61 | } |
| 62 | |
| 63 | void writeNRF(char const* filename, std::vector<SRFPointSource> const& sources, Map const& map, bool normalizeOnset) |
| 64 | { |
| 65 | unsigned numSources = sources.size(); |
| 66 | double minTinit = std::numeric_limits<double>::max(); |
| 67 | unsigned numSamples[3] = { 0, 0, 0 }; |
| 68 | |
| 69 | for (std::vector<SRFPointSource>::const_iterator source = sources.begin(); source != sources.end(); ++source) { |
| 70 | minTinit = std::min(minTinit, source->tinit); |
| 71 | for (unsigned sr = 0; sr < 3; ++sr) { |
| 72 | numSamples[sr] += source->slipRate[sr].size(); |
| 73 | } |
| 74 | } |
| 75 | if (normalizeOnset) { |
| 76 | std::cout << "Minimal tinit: " << minTinit << " -> 0" << std::endl; |
| 77 | } else { |
| 78 | minTinit = 0.0; |
| 79 | } |
| 80 | |
| 81 | Offsets* offsets = new Offsets[numSources+1]; |
| 82 | memset(offsets, 0, (numSources+1) * sizeof(Offsets)); |
| 83 | |
| 84 | double** sliprates = new double*[3]; |
| 85 | for (unsigned sr = 0; sr < 3; ++sr) { |
| 86 | sliprates[sr] = new double[numSamples[sr]]; |
| 87 | memset(sliprates[sr], 0, numSamples[sr] * sizeof(double)); |
| 88 | } |
| 89 | |
| 90 | Point3d* centres = new Point3d[numSources]; |
| 91 | Subfault* subfaults = new Subfault[numSources]; |
| 92 | for (unsigned i = 0; i < numSources; ++i) { |
| 93 | SRFPointSource const& source = sources[i]; |
| 94 | Point3d& centre = centres[i]; |
| 95 | Subfault& sf = subfaults[i]; |
| 96 | sf.tinit = source.tinit - minTinit; |
| 97 | sf.timestep = source.dt; |
| 98 | |
| 99 | #ifdef noproj |
| 100 | centre.x = source.longitude; |
| 101 | centre.y = source.latitude; |
| 102 | centre.z = source.depth; |
| 103 | #else |
| 104 | map.map(source.longitude, source.latitude, source.depth, ¢re.x, ¢re.y, ¢re.z); |
| 105 | #endif |
| 106 | |
| 107 | map.toMCS(source.strike, source.dip, source.rake, 1.0, 0.0, 0.0, &sf.tan1(0), &sf.tan1(1), &sf.tan1(2)); |
| 108 | map.toMCS(source.strike, source.dip, source.rake, 0.0, 1.0, 0.0, &sf.tan2(0), &sf.tan2(1), &sf.tan2(2)); |
| 109 | map.toMCS(source.strike, source.dip, source.rake, 0.0, 0.0, 1.0, &sf.normal(0), &sf.normal(1), &sf.normal(2)); |
| 110 | |
| 111 | // g / (cm s^2) -> Pa (= kg / (m s^2)) |
| 112 | sf.mu = source.shearModulus * 1.0e-1; |
| 113 | // cm^2 -> m^2 |
| 114 | sf.area = source.area * 1.0e-4; |
| 115 | for (unsigned sr = 0; sr < 3; ++sr) { |
| 116 | unsigned offset = offsets[i][sr]; |
| 117 | std::copy(source.slipRate[sr].begin(), source.slipRate[sr].end(), &sliprates[sr][ offset ]); |
| 118 | offsets[i+1][sr] = offset + source.slipRate[sr].size(); |
| 119 | } |
| 120 | } |