00001 /*************************************************************************** 00002 ** 00003 ** This file is part of QGpCoreWave. 00004 ** 00005 ** This library is free software; you can redistribute it and/or 00006 ** modify it under the terms of the GNU Lesser General Public 00007 ** License as published by the Free Software Foundation; either 00008 ** version 2.1 of the License, or (at your option) any later version. 00009 ** 00010 ** This file is distributed in the hope that it will be useful, but WITHOUT 00011 ** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 ** FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 00013 ** License for more details. 00014 ** 00015 ** You should have received a copy of the GNU Lesser General Public 00016 ** License along with this library; if not, write to the Free Software 00017 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00018 ** 00019 ** See http://www.geopsy.org for more information. 00020 ** 00021 ** Created : 2004-10-21 00022 ** Authors : 00023 ** Marc Wathelet 00024 ** Marc Wathelet (ULg, Liège, Belgium) 00025 ** Marc Wathelet (LGIT, Grenoble, France) 00026 ** 00027 ***************************************************************************/ 00028 00029 #ifndef DISPERSION_H 00030 #define DISPERSION_H 00031 00032 #include <math.h> 00033 00034 #include <QGpCoreTools.h> 00035 00036 #include "Ellipticity.h" 00037 #include "ModalStorage.h" 00038 #include "Love.h" 00039 #include "Rayleigh.h" 00040 #include "QGpCoreWaveDLLExport.h" 00041 00042 namespace QGpCoreWave { 00043 00044 class Dispersion : public ModalStorage 00045 { 00046 TRANSLATIONS("Dispersion") 00047 public: 00048 QGPCOREWAVE_EXPORT Dispersion(int nModes, const QVector<double>* x); 00049 00050 QGPCOREWAVE_EXPORT void setGroupSlowness(); 00051 00052 bool calculate(Rayleigh * waveModel, Ellipticity * ell=0, QAtomicInt * terminated=0) 00053 { 00054 return calculate<double, Rayleigh>(waveModel, ell, terminated); 00055 } 00056 #ifdef BIGNUM 00057 bool calculate(RayleighBigNum * waveModel, Ellipticity * ell=0, QAtomicInt * terminated=0) 00058 { 00059 return calculate<mp_real, RayleighBigNum>(waveModel, ell, terminated); 00060 } 00061 #endif 00062 bool calculate(Love * waveModel, QAtomicInt * terminated=0) 00063 { 00064 return calculate<double, Love>(waveModel, 0, terminated); 00065 } 00066 void setPrecision(double p) {_precision=p;} 00067 QGPCOREWAVE_EXPORT bool refine(int modeIndex, int lastIndex, double omega, RootSolver<Rayleigh>& modelS, Ellipticity& ell); 00068 QGPCOREWAVE_EXPORT ModalStorage * deltaK() const; 00069 private: 00070 template <class RealType, class WaveModel> 00071 bool calculate(WaveModel * waveModel, Ellipticity * ell, QAtomicInt * terminate); 00072 double _precision; 00073 }; 00074 00075 template <class RealType, class WaveModel> 00076 bool Dispersion::calculate(WaveModel * waveModel, Ellipticity * ell, QAtomicInt * terminated) 00077 { 00078 TRACE; 00079 RealType * curMode=new RealType [ xCount() ]; 00080 RealType * lastMode=new RealType [ xCount() ]; 00081 if(ell) ell->refineClear(); 00082 refineClear(); 00083 //#define calculate_DEBUG 00084 00085 const Seismic1DModel * model=waveModel->model(); 00086 RootSolverTemplate<RealType, WaveModel> solver(waveModel); 00087 RealType minSlow=model->minSlowS(); 00088 RealType maxSlow, curSlow, lastSupBound=model->maxSlowR(), newCurSlow; 00089 bool forbitVelocityInversion=model->checkVelocityInversion()==-1; 00090 // Initialize equation's polarity 00091 waveModel->setOmega(0.05); // the sign of polarity is the same at all frequencies 00092 solver.setPolarity(model->maxSlowR()); 00093 int lastOmega=xCount() - 1; 00094 00095 for(int iMode=0; iMode < modeCount(); iMode++ ) { 00096 #ifdef calculate_DEBUG 00097 printf( "***************** Mode %i ******************\n", iMode); 00098 #endif 00099 if(ell) { 00100 solver.setPrecision(_precision*1e-3); // Down to 1e-10 in default case (see constructor) 00101 } else { 00102 solver.setPrecision(_precision); // Down to 1e-7 in default case (see constructor) 00103 } 00104 if(iMode==0) { 00105 maxSlow=model->maxSlowR(); 00106 // 2009-02-15 Reduce factor from 0.5 to 0.25 to avoid too many tries 00107 // With 0.5 it appears that mode jumping occurs quite a lot 00108 solver.setRelativeStep(0.25 * (model->maxSlowR() - model->maxSlowS())/maxSlow); 00109 } else { 00110 RealType * tmpMode=lastMode; 00111 lastMode=curMode; 00112 curMode=tmpMode; 00113 maxSlow=lastMode[ lastOmega ]; 00114 solver.setRelativeStep(0.5 * (model->maxSlowR() - model->maxSlowS())/maxSlow * (2 * M_PI/x(lastOmega) )); 00115 solver.inversePolarity(); 00116 } 00117 RealValue * ellValues; 00118 if(ell && iMode<ell->modeCount()) { 00119 ellValues=ell->mode(iMode); 00120 } else { 00121 ellValues=0; 00122 } 00123 RealValue * slownessValues=mode(iMode); 00124 00125 while(true) { 00126 bool errorYetDetected=false; 00127 curSlow=maxSlow; 00128 for(int i=lastOmega;i >= 0;i-- ) { 00129 #ifdef calculate_DEBUG 00130 printf( "***************** Frequency %i omega %lf f %lf\n", i, x(i), x(i)/(2*M_PI)); 00131 #endif 00132 if(terminated && terminated->testAndSetOrdered(true,true)) { 00133 // Catch user interrupt 00134 delete [] curMode; 00135 delete [] lastMode; 00136 return false; 00137 } 00138 waveModel->setOmega(x( i) ); 00139 if(iMode > 0) { 00140 maxSlow=lastMode[ i ]; 00141 if(curSlow > maxSlow) curSlow=maxSlow; 00142 } 00143 //printf("%20.10lf %20.10lf %20.10lf\n",dble(curSlow),dble(_minSlow),dble(maxSlow)); 00144 if(solver.searchDown(curSlow, minSlow, maxSlow) ) { 00145 solver.neville(); 00146 newCurSlow=solver.lower(); 00147 } else { 00148 if(iMode==0) { 00149 // for fundamental mode this condition can never occur unless a mode jumping occured 00150 errorYetDetected=true; 00151 App::stream() << tr( "** Warning ** : mode jumping for mode %1 (end), reducing " 00152 "step ratio to %2" ).arg(iMode).arg(Number::toDouble(solver.searchStep() * 0.1) ) << endl; 00153 break; 00154 } else { 00155 // The end of the curve, the mode does not exist at this period (no root) 00156 for(int j=i;j >= 0;j-- ) { 00157 slownessValues[ j ].setValid(false); 00158 curMode[ j ]=minSlow; 00159 } 00160 if(ellValues) { 00161 for(int j=i;j >= 0;j-- ) ellValues[ j ].setValid(false); 00162 } 00163 curSlow=minSlow; 00164 break; 00165 } 00166 } 00167 if(ellValues) { 00168 /* Take a mid slowness between the refined bounds and recompute the 00169 propagator matrices to get the corresponding sub-determinents necessary 00170 for ellipticity computation. Keep ellipticity with its sign. */ 00171 waveModel->y(0.5*(solver.upper()+newCurSlow)); 00172 ellValues[ i ].setValue(waveModel->ellipticity()); 00173 } 00174 double newLastSupBound=solver.upper(); // searchUp may alter it 00175 //printf("%20.10lf %20.10lf %20.10lf\n",x(i)/(2*M_PI), modelS.infBoundRoot(), modelS.supBoundRoot()); 00176 if(i < lastOmega && newCurSlow - curSlow > _precision) { 00177 if(forbitVelocityInversion) { 00178 errorYetDetected=true; 00179 App::stream() << tr( "** Warning ** : retrograde dispersion not allowed (mode %1), " 00180 "reducing step ratio to %2" ).arg(iMode).arg(Number::toDouble(solver.searchStep() * 0.1) ) << endl; 00181 break; 00182 } else { 00183 // Inversions are allowed, check at last point that no mode exist at a higher slowness 00184 waveModel->setOmega(x( i + 1) ); 00185 solver.setRelativeStep(solver.searchStep() * 0.1); 00186 if(solver.searchUp(lastSupBound, iMode>0 ? lastMode[ i + 1 ] : maxSlow) ) { 00187 // recalculating curve because of mode jumping 00188 App::stream() << tr( "** Warning ** : mode jumping for mode %1 (middle at %3 Hz), " 00189 "reducing step ratio to %2" ) 00190 .arg(iMode) 00191 .arg(Number::toDouble(solver.searchStep()) ) 00192 .arg(x( i)/ (2*M_PI)) << endl; 00193 errorYetDetected=true; 00194 solver.setRelativeStep(solver.searchStep() * 10.0); 00195 break; 00196 } 00197 solver.setRelativeStep(solver.searchStep() * 10.0); 00198 } 00199 } 00200 lastSupBound=newLastSupBound; 00201 curSlow=newCurSlow; 00202 curMode[ i ]=curSlow; 00203 slownessValues[ i ].setValue(Number::toDouble(curSlow) ); 00204 slownessValues[ i ].setValid(true); 00205 #ifdef calculate_DEBUG 00206 printf( "%30.20lf %30.20lf\n", 0.5 * x(i)/M_PI, slownessValues[ i ].value()); 00207 #endif 00208 } 00209 // When calculating the dispersion curve, stepRatio is positive to search downwards 00210 // (From maxSlow towards minSlow) 00211 // Here we want to search carefully upwards: from curSlow to maxSlow 00212 // If at least one root is found, this proves that mode jumping occured 00213 if(curSlow > minSlow) { 00214 // The curSlow was calculated using _solve and _x2 is the value returned by 00215 // Solve (slightly greater than the true root). 00216 // Here we need the _x1 value, e.g. slightly less than the true root, in 00217 // order to be sure that the root eventually found will not be the same 00218 // as curSlow 00219 curSlow=solver.upper(); 00220 } 00221 if(errorYetDetected || solver.searchUp(curSlow, maxSlow) ) { 00222 // recalculating curve because of mode jumping 00223 if(solver.searchStep() < 1e-3) { 00224 App::stream() << tr( "** Warning ** : dispersion curve impossible to obtain " 00225 "for mode %1" ).arg(iMode) << endl; 00226 for( ; iMode < modeCount(); iMode++ ) { 00227 RealValue * slownessValues=mode(iMode); 00228 for(int j=lastOmega;j >= 0;j-- ) { 00229 slownessValues[ j ].setValid(false); 00230 } 00231 } 00232 delete [] curMode; 00233 delete [] lastMode; 00234 return false; 00235 } 00236 solver.setPrecision(solver.precision() * 0.1); 00237 solver.setRelativeStep(solver.searchStep() * 0.1); 00238 if( !errorYetDetected) 00239 App::stream() << tr( "** Warning ** : mode jumping for mode %1 (end), reducing " 00240 "step ratio to %2" ).arg(iMode).arg(Number::toDouble(solver.searchStep()) ) << endl; 00241 if(iMode > 0) maxSlow=lastMode[ lastOmega ]; 00242 } else break; 00243 } 00244 } 00245 delete [] curMode; 00246 delete [] lastMode; 00247 return true; 00248 } 00249 00250 } // namespace QGpCoreWave 00251 00252 #endif // DISPERSION_H
