| 1 | #include "machdefs.h"
|
|---|
| 2 | #include <stdio.h>
|
|---|
| 3 | #include <stdlib.h>
|
|---|
| 4 | #include <math.h>
|
|---|
| 5 |
|
|---|
| 6 | #include "fsvcache.h"
|
|---|
| 7 | #include "nbtri.h"
|
|---|
| 8 | #include "nbrandom.h"
|
|---|
| 9 | #include "nbmath.h"
|
|---|
| 10 | #include "nbinteg.h"
|
|---|
| 11 | #include "strutil.h"
|
|---|
| 12 | #include "fsvst.h"
|
|---|
| 13 | #include "nbsread.h"
|
|---|
| 14 | #include "nbgene.h"
|
|---|
| 15 |
|
|---|
| 16 | /*==========================================================================*/
|
|---|
| 17 | int Give_Annee_Obs(int lp)
|
|---|
| 18 | /* retourne l'annee d'observation codee en entier du fichier de suivi utilise */
|
|---|
| 19 | {
|
|---|
| 20 | int j,anneeobs;
|
|---|
| 21 | double t=0.;
|
|---|
| 22 | anneeobs = -1;
|
|---|
| 23 | for(j=0;j<nmes[0];j++) {
|
|---|
| 24 | t = timeu[0][j].TStart / JourSec;
|
|---|
| 25 | if( t <= 0. ) continue;
|
|---|
| 26 | /* confusion possible smc 93 avec debut lmc 9394 */
|
|---|
| 27 | if( t > 1288. && t < 1314. ) continue;
|
|---|
| 28 | if( t > 700. && t < 850. ) {anneeobs=9192; break;}
|
|---|
| 29 | if( t > 950. && t < 1200. ) {anneeobs=9293; break;}
|
|---|
| 30 | if( t > 1300. && t < 1570. ) {anneeobs=9394; break;}
|
|---|
| 31 | if( t > 1630. && t < 1900. ) {anneeobs=9495; break;}
|
|---|
| 32 | if( t > 2300. && t < 9999. ) {anneeobs=9697; break;}
|
|---|
| 33 | }
|
|---|
| 34 | if(lp>0) printf("annee_obs = %d (j=%d t=%f)\n",anneeobs,j,t);
|
|---|
| 35 | return (anneeobs);
|
|---|
| 36 | }
|
|---|
| 37 |
|
|---|
| 38 | /*=========================================================================*/
|
|---|
| 39 | float Rayon_SP_Th(float magv)
|
|---|
| 40 | /* Fournit le rayon de l'etoile de SP en rayons solaires
|
|---|
| 41 | * en fonction de sa magnitude Vj.
|
|---|
| 42 | * Formule theorique uniquement pour la sequence principale:
|
|---|
| 43 | * Rstar_SP = exp(-0.22 * (magv - MagV_sol) ) avec MagV_sol = 5.1
|
|---|
| 44 | */
|
|---|
| 45 | {
|
|---|
| 46 | return( (float) exp(-0.22*(magv-5.1)) );
|
|---|
| 47 | }
|
|---|
| 48 |
|
|---|
| 49 | /*==========================================================================*/
|
|---|
| 50 | void AMPLML(int io)
|
|---|
| 51 | /* Generation MonteCarlo */
|
|---|
| 52 | {
|
|---|
| 53 | int i,ic,j,tf;
|
|---|
| 54 | float x;
|
|---|
| 55 | double t,tt,u2,u;
|
|---|
| 56 | int Nseq_tf = 3;
|
|---|
| 57 | unsigned short seed_16v[3];
|
|---|
| 58 |
|
|---|
| 59 | mc.TauSim = mc.TauSimI;
|
|---|
| 60 | mc.T0Sim = mc.T0SimI;
|
|---|
| 61 | mc.U0Sim = mc.U0SimI;
|
|---|
| 62 | mc.USim = mc.USimI;
|
|---|
| 63 |
|
|---|
| 64 | if(io == 1) {
|
|---|
| 65 |
|
|---|
| 66 | /* initialisation des generateurs */
|
|---|
| 67 | seed_16v[0]=(unsigned short) mc.iseed1;
|
|---|
| 68 | seed_16v[1]=(unsigned short) mc.iseed2;
|
|---|
| 69 | seed_16v[2]=(unsigned short) mc.iseed3;
|
|---|
| 70 | Ini_Ranf(seed_16v,0);
|
|---|
| 71 |
|
|---|
| 72 | mc.NGene = mc.NPoints_ampli = mc.NPoints_ampli_tf = mc.Net_ampli_tf = 0;
|
|---|
| 73 |
|
|---|
| 74 | if(mc.montecar<=0) return;
|
|---|
| 75 | if ( mc.Taille_Finie > 2 )
|
|---|
| 76 | { printf("mauvaise valeur mc.Taille_Finie %d\n",mc.Taille_Finie);
|
|---|
| 77 | exit(-1); }
|
|---|
| 78 | if ( mc.Blending > 2 || mc.Blending < 0 )
|
|---|
| 79 | { printf("mauvaise valeur mc.Blending %d\n",mc.Blending);
|
|---|
| 80 | exit(-1); }
|
|---|
| 81 |
|
|---|
| 82 | printf("MONTECARLO %2d: seed %d %d %d\n",io,mc.iseed1,mc.iseed2,mc.iseed3);
|
|---|
| 83 | printf(" Generation %d, Taille_Finie %d blending %d\n"
|
|---|
| 84 | ,mc.montecar,mc.Taille_Finie,mc.Blending);
|
|---|
| 85 | printf(" Nombre de generations par etoile %d\n",mc.NGenStar);
|
|---|
| 86 | printf(" TauSim=%f TauMinI=%f TauMaxI=%f TireTau=%p\n"
|
|---|
| 87 | ,mc.TauSimI,mc.TauMinI,mc.TauMaxI,mc.TireTau);
|
|---|
| 88 | printf(" T0Sim=%f T0MinI=%f T0MaxI=%f TireT0=%p\n"
|
|---|
| 89 | ,mc.T0SimI,mc.T0MinI,mc.T0MaxI,mc.TireT0);
|
|---|
| 90 | printf(" U0Sim=%f U0MinI=%f U0MaxI=%f TireU0=%p\n"
|
|---|
| 91 | ,mc.U0SimI,mc.U0MinI,mc.U0MaxI,mc.TireU0);
|
|---|
| 92 | printf(" USim =%f UMinI =%f UMaxI =%f TireU=%p\n"
|
|---|
| 93 | ,mc.USimI,mc.UMinI,mc.UMaxI,mc.TireU);
|
|---|
| 94 | printf(" coupure pour taille finie a U/u0 < %f\n",mc.Usu0_Cut);
|
|---|
| 95 |
|
|---|
| 96 | if( (mc.U0MaxI <= mc.U0MinI && mc.U0SimI <= 0. && mc.Taille_Finie==0 )
|
|---|
| 97 | || (mc.U0MaxI <= mc.U0MinI && mc.U0SimI < 0. && mc.Taille_Finie>0 ) ) {
|
|---|
| 98 | printf("Mauvaises valeurs U0SimI U0MinI U0MaxI\n");
|
|---|
| 99 | exit(-1);
|
|---|
| 100 | }
|
|---|
| 101 | if( mc.UMaxI <= mc.UMinI && mc.USimI <= 0. && mc.Taille_Finie>0 ) {
|
|---|
| 102 | printf("Mauvaises valeurs USimI UMinI UMaxI\n");
|
|---|
| 103 | exit(-1);
|
|---|
| 104 | }
|
|---|
| 105 | if( mc.TauMaxI <= mc.TauMinI && mc.TauSimI <= 0.) {
|
|---|
| 106 | printf("Mauvaises valeurs TauSimI TauMinI TauMaxI\n");
|
|---|
| 107 | exit(-1);
|
|---|
| 108 | }
|
|---|
| 109 |
|
|---|
| 110 | } else if(io == 2) {
|
|---|
| 111 |
|
|---|
| 112 | if( mc.T0MaxI <= mc.T0MinI && mc.T0SimI <= 0. ) {
|
|---|
| 113 | mc.T0MinI_eff=TFIRST;
|
|---|
| 114 | mc.T0MaxI_eff=TLAST;
|
|---|
| 115 | } else {
|
|---|
| 116 | mc.T0MinI_eff=mc.T0MinI;
|
|---|
| 117 | mc.T0MaxI_eff=mc.T0MaxI;
|
|---|
| 118 | }
|
|---|
| 119 |
|
|---|
| 120 | if(mc.montecar<=0) return;
|
|---|
| 121 |
|
|---|
| 122 | if( mc.Blending>0 && mc.Get_Par_Blending==0 ) {
|
|---|
| 123 | printf("Vous avez demande une generation avec du Blending\n");
|
|---|
| 124 | printf(" mais vous n'avez pas fourni la generation des parametres\n");
|
|---|
| 125 | exit(-1);
|
|---|
| 126 | }
|
|---|
| 127 |
|
|---|
| 128 | printf("MONTECARLO %2d: ccd %6d, T0MinI_eff=%f T0MaxI_eff=%f\n"
|
|---|
| 129 | ,io,ccdnum,mc.T0MinI_eff,mc.T0MaxI_eff);
|
|---|
| 130 |
|
|---|
| 131 | } else if(io == 3) {
|
|---|
| 132 |
|
|---|
| 133 | mc.NGene++;
|
|---|
| 134 |
|
|---|
| 135 | /* Blending preparation */
|
|---|
| 136 | if( mc.Blending>0 ) mc.Get_Par_Blending();
|
|---|
| 137 |
|
|---|
| 138 | /* Tirage aleatoire de u0 */
|
|---|
| 139 | if ( mc.TireU0 != 0 ) {
|
|---|
| 140 | x = drand01();
|
|---|
| 141 | mc.U0Sim = mc.TireU0(x);
|
|---|
| 142 | } else if(mc.U0MaxI > mc.U0MinI) {
|
|---|
| 143 | mc.U0Sim= -1.;
|
|---|
| 144 | while( mc.U0Sim <= 0. ) {
|
|---|
| 145 | mc.U0Sim = mc.U0MinI + (mc.U0MaxI-mc.U0MinI)*drand01();
|
|---|
| 146 | if(mc.U0Sim<=0.) printf("ATTENTION: U0 SIMULE A %f min=%f max=%f\n"
|
|---|
| 147 | ,mc.U0Sim,mc.U0MinI,mc.U0MaxI);
|
|---|
| 148 | } }
|
|---|
| 149 |
|
|---|
| 150 | /* generation plate de tau */
|
|---|
| 151 | if ( mc.TireTau != 0 ) {
|
|---|
| 152 | x = drand01();
|
|---|
| 153 | mc.TauSim = mc.TireTau(x);
|
|---|
| 154 | } else if ( mc.TauMaxI > mc.TauMinI ) {
|
|---|
| 155 | mc.TauSim = mc.TauMinI
|
|---|
| 156 | + (mc.TauMaxI-mc.TauMinI)*drand01();
|
|---|
| 157 | }
|
|---|
| 158 |
|
|---|
| 159 | /* Tirage aleatoire de t0 */
|
|---|
| 160 | if ( mc.TireT0 != 0 ) {
|
|---|
| 161 | x = drand01();
|
|---|
| 162 | mc.T0Sim = mc.TireT0(x);
|
|---|
| 163 | } else if( mc.T0MaxI_eff > mc.T0MinI_eff ) {
|
|---|
| 164 | mc.T0Sim = mc.T0MinI_eff
|
|---|
| 165 | + (mc.T0MaxI_eff-mc.T0MinI_eff)*drand01();
|
|---|
| 166 | }
|
|---|
| 167 |
|
|---|
| 168 | /* generation plate de U */
|
|---|
| 169 | if( mc.Taille_Finie > 0 ) {
|
|---|
| 170 | if ( mc.TireU != 0 ) {
|
|---|
| 171 | x = drand01();
|
|---|
| 172 | mc.USim = mc.TireU(x);
|
|---|
| 173 | } else if ( mc.Taille_Finie>0 && mc.UMaxI > mc.UMinI ) {
|
|---|
| 174 | mc.USim= mc.UMinI + (mc.UMaxI-mc.UMinI)*drand01();
|
|---|
| 175 | }
|
|---|
| 176 | if(mc.Taille_Finie==2) {
|
|---|
| 177 | if( mc.Blending==0 ) mc.USim *= starcal.Rstar;
|
|---|
| 178 | else if( mc.Blending==1 ) mc.USim *= mc.starcal1.Rstar;
|
|---|
| 179 | else if( mc.Blending==2 ) mc.USim *= mc.starcal2.Rstar;
|
|---|
| 180 | }
|
|---|
| 181 | }
|
|---|
| 182 |
|
|---|
| 183 | /* calcul amplification ponctuelle/taille finie (paczinsky) */
|
|---|
| 184 | mc.A0MaxPt = GRAND;
|
|---|
| 185 | if( mc.U0Sim > 0. ) mc.A0MaxPt = ampli_ponct(mc.U0Sim);
|
|---|
| 186 | mc.A0Max = (mc.Taille_Finie && mc.USim>0.) ?
|
|---|
| 187 | ampli_tf(mc.U0Sim,mc.USim,Nseq_tf) : mc.A0MaxPt;
|
|---|
| 188 | if( mc.Blending>0) for(ic=0;ic<NCOULMX;ic++) {
|
|---|
| 189 | if( mc.Blending == 1 )
|
|---|
| 190 | mc.A0Max_blend[ic] = (mc.A0Max -1.)* mc.coeff_Arec1[ic] + 1.;
|
|---|
| 191 | else if( mc.Blending == 2 )
|
|---|
| 192 | mc.A0Max_blend[ic] = (mc.A0Max -1.)* mc.coeff_Arec2[ic] + 1.;
|
|---|
| 193 | }
|
|---|
| 194 |
|
|---|
| 195 | tf = 0;
|
|---|
| 196 | for(ic=0;ic<NCOULMX;ic++) {
|
|---|
| 197 | for(i=0;i<nmes[ic];i++) {
|
|---|
| 198 | t = date[ic][i];
|
|---|
| 199 | tt = (t-mc.T0Sim)/mc.TauSim;
|
|---|
| 200 | u2 = mc.U0Sim*mc.U0Sim + tt*tt;
|
|---|
| 201 | u = sqrt(u2);
|
|---|
| 202 | ampli[ic][i] = ampli_ponct(u);
|
|---|
| 203 | mc.NPoints_ampli++;
|
|---|
| 204 | if( mc.Taille_Finie && t<GRAND2 && mc.USim>0. ) {
|
|---|
| 205 | if( mc.USim > mc.Usu0_Cut * u ) {
|
|---|
| 206 | ampli[ic][i] = ampli_tf(u,mc.USim,Nseq_tf);
|
|---|
| 207 | mc.NPoints_ampli_tf++;
|
|---|
| 208 | if(tf==0) { mc.Net_ampli_tf++ ; tf = 1;}
|
|---|
| 209 | }
|
|---|
| 210 | }
|
|---|
| 211 | if(mc.Blending==1)
|
|---|
| 212 | ampli[ic][i] = (ampli[ic][i] -1.)* mc.coeff_Arec1[ic] + 1.;
|
|---|
| 213 | else if(mc.Blending==2)
|
|---|
| 214 | ampli[ic][i] = (ampli[ic][i] -1.)* mc.coeff_Arec2[ic] + 1.;
|
|---|
| 215 | } }
|
|---|
| 216 |
|
|---|
| 217 | if(debug>1) {
|
|---|
| 218 | for(ic=0;ic<NCOULMX;ic++) {
|
|---|
| 219 | printf("ampli[%3d]-1 u0=%12.5g t0=%12.5g tau=%12.5g U=%12.5g:\n"
|
|---|
| 220 | ,ic,mc.U0Sim,mc.T0Sim,mc.TauSim,mc.USim);
|
|---|
| 221 | for(i=0;i<nmesure[ic];i++)
|
|---|
| 222 | { j=indexu[ic][i]; printf(" %12.5g",ampli[ic][j]-1.);
|
|---|
| 223 | if(i%10==9) printf("\n");}
|
|---|
| 224 | if((nmesure[ic]-1)%10!=9) printf("\n");
|
|---|
| 225 | } }
|
|---|
| 226 |
|
|---|
| 227 | } else if(io == 4) {
|
|---|
| 228 |
|
|---|
| 229 | return;
|
|---|
| 230 |
|
|---|
| 231 | } else if(io == 5) {
|
|---|
| 232 |
|
|---|
| 233 | if(mc.montecar<=0) return;
|
|---|
| 234 |
|
|---|
| 235 | printf("MONTECARLO %2d: Nombre d'evenements engendres %d\n",io,mc.NGene);
|
|---|
| 236 | printf("Nb de pts amplifies %d, nb de pts amplifies avec taille finie %d\n"
|
|---|
| 237 | ,mc.NPoints_ampli,mc.NPoints_ampli_tf);
|
|---|
| 238 | printf("Nb d'etoiles amplifiees avec taille finie %d\n"
|
|---|
| 239 | ,mc.Net_ampli_tf);
|
|---|
| 240 | Get_Ranf(seed_16v,0);
|
|---|
| 241 | mc.iseed1 = (int_4)seed_16v[0];
|
|---|
| 242 | mc.iseed2 = (int_4)seed_16v[1];
|
|---|
| 243 | mc.iseed3 = (int_4)seed_16v[2];
|
|---|
| 244 | printf("Aleatoire_de_fin_de_Generation %d %d %d\n"
|
|---|
| 245 | ,mc.iseed1,mc.iseed2,mc.iseed3);
|
|---|
| 246 |
|
|---|
| 247 | }
|
|---|
| 248 | return;
|
|---|
| 249 | }
|
|---|
| 250 |
|
|---|
| 251 | /*==========================================================================*/
|
|---|
| 252 | double ampli_ponct(double u)
|
|---|
| 253 | /* amplification paczinsky dans le cas d'une source ponctuelle */
|
|---|
| 254 | {
|
|---|
| 255 | u *= u;
|
|---|
| 256 | if(u==0.) return(GRAND2);
|
|---|
| 257 | u = (u+2.)/sqrt(u*(u+4.));
|
|---|
| 258 | return (u);
|
|---|
| 259 | }
|
|---|
| 260 |
|
|---|
| 261 | /*==========================================================================*/
|
|---|
| 262 | double ampli_tf(double u, double U,int Nseq)
|
|---|
| 263 | /*
|
|---|
| 264 | Amplification paczinsky dans le cas d'une source non ponctuelle.
|
|---|
| 265 | ( integration analytique partie interieure et
|
|---|
| 266 | integration numerique partie exterieure )
|
|---|
| 267 | Integration exterieure par methode de Gauss (Higher order terms).
|
|---|
| 268 | u = parametre d'impact en rayon d'Einstein en fct du temps
|
|---|
| 269 | U = rayon de l'etoile projete en rayon d'Einstein
|
|---|
| 270 | Nseq = nombre de pas pour l'integrale
|
|---|
| 271 | */
|
|---|
| 272 | {
|
|---|
| 273 | int i,j;
|
|---|
| 274 | double u2,U2,lim1,lim2,dlim,A0,xc,x,nx,b,teta,Ka,z;
|
|---|
| 275 | int mIhoqN;
|
|---|
| 276 | double *IhoqNX, *IhoqNW;
|
|---|
| 277 |
|
|---|
| 278 | if( Nseq <= 0 ) Nseq = 2;
|
|---|
| 279 | if( u<0. || U<0. ) return(-1.);
|
|---|
| 280 | if( u==0. && U==0. ) return(GRAND2);
|
|---|
| 281 |
|
|---|
| 282 | mIhoqN=mIhoq9; IhoqNX=Ihoq9X; IhoqNW=Ihoq9W;
|
|---|
| 283 |
|
|---|
| 284 | u2 = u*u;
|
|---|
| 285 | U2 = U*U;
|
|---|
| 286 |
|
|---|
| 287 | Ka = 0.;
|
|---|
| 288 | lim1 = -1.;
|
|---|
| 289 | lim2 = 1.;
|
|---|
| 290 | dlim = (lim2 - lim1)/Nseq;
|
|---|
| 291 |
|
|---|
| 292 | /*** CAS 1 : u<=U (u=0 possible si U!=0) *** */
|
|---|
| 293 | /* A = A0 + A1 */
|
|---|
| 294 | /* A0 = partie integrable (sum(dx,0<x<U-u)) */
|
|---|
| 295 | /* = (U-u) * sqrt((U-u)**2 + 4) / U2 */
|
|---|
| 296 | /* A1 = 2*u/(Pi*U2) * sum(dx,-1<x<1) */
|
|---|
| 297 | /* {(u2+2)/sqrt(u2+4)*acos[(u*(1+x**2)/2+U*x)/(U+x*u)]} */
|
|---|
| 298 |
|
|---|
| 299 | if(u<=U) {
|
|---|
| 300 | /* on calcule d'abord la partie analytique A0 */
|
|---|
| 301 | x = (U-u);
|
|---|
| 302 | A0 = x * sqrt(x*x + 4.);
|
|---|
| 303 |
|
|---|
| 304 | /* on integre ensuite la partie numerique A1 */
|
|---|
| 305 | for(j=0;j<Nseq;j++) {
|
|---|
| 306 | xc = lim1 + (j+0.5)*dlim;
|
|---|
| 307 | for(i=0;i<mIhoqN;i++) {
|
|---|
| 308 | x = xc + IhoqNX[i]*dlim;
|
|---|
| 309 | nx = U + x*u;
|
|---|
| 310 | teta = acos( ( u*(1.+x*x)/2. + U*x ) /nx );
|
|---|
| 311 | nx *= nx;
|
|---|
| 312 | b = (nx + 2.) / sqrt(nx + 4.);
|
|---|
| 313 | Ka += b * teta * IhoqNW[i];
|
|---|
| 314 | } }
|
|---|
| 315 | A0 += 2. * u / Pi * Ka * dlim;
|
|---|
| 316 | A0 /= U2;
|
|---|
| 317 | }
|
|---|
| 318 |
|
|---|
| 319 | /*** CAS 2 : u>U --> U peut tendre vers 0 *** */
|
|---|
| 320 | /* A = 2/(Pi*U) * sum(dx,-1<x<1) */
|
|---|
| 321 | /* {(u2+2)/sqrt(u2+4)*acos[1-U2*(1-x**2)/(2*u*(u+x*U))]} */
|
|---|
| 322 | /* quand U->0, la divergence (1/U) est enlevee en developpant */
|
|---|
| 323 | /* acos(1-z) = sqrt(2*z)*[1+z/12+3*z2/160+...] */
|
|---|
| 324 |
|
|---|
| 325 | else if(u>U) {
|
|---|
| 326 | for(j=0;j<Nseq;j++) {
|
|---|
| 327 | xc = lim1 + (j+0.5)*dlim;
|
|---|
| 328 | for(i=0;i<mIhoqN;i++) {
|
|---|
| 329 | x = xc + IhoqNX[i]*dlim;
|
|---|
| 330 | nx = u + x*U;
|
|---|
| 331 | z = (1.-x*x) /(2.*u*nx);
|
|---|
| 332 | nx *= nx;
|
|---|
| 333 | b = (nx + 2.) / sqrt(nx + 4.);
|
|---|
| 334 | if(U > 0.001*u) {
|
|---|
| 335 | teta = acos( 1. - U2*z ) / U;
|
|---|
| 336 | } else { /* cas U-> 0 */
|
|---|
| 337 | teta = sqrt(2.*z);
|
|---|
| 338 | z *= U2;
|
|---|
| 339 | teta *= 1. + z*(0.08333333333333 + 0.01875*z);
|
|---|
| 340 | }
|
|---|
| 341 | Ka += b * teta * IhoqNW[i];
|
|---|
| 342 | } }
|
|---|
| 343 | A0 = 2. / Pi * Ka * dlim;
|
|---|
| 344 | }
|
|---|
| 345 |
|
|---|
| 346 | /*** CAS IMPOSSIBLE ****/
|
|---|
| 347 | else {
|
|---|
| 348 | printf("ampli_tf: CAS IMPOSSIBLE u=%f U=%f\n",u,U);
|
|---|
| 349 | A0 = -1.;
|
|---|
| 350 | }
|
|---|
| 351 |
|
|---|
| 352 | return (A0);
|
|---|
| 353 | }
|
|---|
| 354 |
|
|---|
| 355 | /*==========================================================================*/
|
|---|
| 356 | double u0_ponct(double A)
|
|---|
| 357 | /*
|
|---|
| 358 | Parametre d'impact pour une amplification A
|
|---|
| 359 | dans le cas d'une source ponctuelle.
|
|---|
| 360 | input: A = amplification
|
|---|
| 361 | */
|
|---|
| 362 | {
|
|---|
| 363 | double u0;
|
|---|
| 364 | if(A<1.) return(-1.);
|
|---|
| 365 | if(A==1.) return(GRAND2);
|
|---|
| 366 | u0 = 2. * ( A/sqrt(A*A-1.) - 1.);
|
|---|
| 367 | u0 = sqrt(u0);
|
|---|
| 368 | return(u0);
|
|---|
| 369 | }
|
|---|
| 370 |
|
|---|
| 371 | /*==========================================================================*/
|
|---|
| 372 | double u0_tf(double A,double U,double prec,int *niter,double *Afinal)
|
|---|
| 373 | /*
|
|---|
| 374 | Parametre d'impact pour une amplification A
|
|---|
| 375 | dans le cas d'une source non ponctuelle.
|
|---|
| 376 | input:
|
|---|
| 377 | A = amplification
|
|---|
| 378 | U = Rstar projete / Reinst
|
|---|
| 379 | prec = precision souhaitee sur u0
|
|---|
| 380 | niter = nombre maximum d'iterations permises
|
|---|
| 381 | output:
|
|---|
| 382 | niter = nombre d'iterations utilisees (-1 si pas converge)
|
|---|
| 383 | Afinal = amplification finale calculee pour u0
|
|---|
| 384 | return = u0 si trouve
|
|---|
| 385 | -1. si A<1.
|
|---|
| 386 | -2. si A=1.
|
|---|
| 387 | -3. si U<=0
|
|---|
| 388 | -4. si A trop grand
|
|---|
| 389 | -5. si encadrement de A non trouve
|
|---|
| 390 | */
|
|---|
| 391 | {
|
|---|
| 392 | int nitermx;
|
|---|
| 393 | double a0=0,a1,a2,u0=0,u1,u2,test;
|
|---|
| 394 |
|
|---|
| 395 | nitermx = *niter;
|
|---|
| 396 | if(nitermx<=0) nitermx = 25;
|
|---|
| 397 |
|
|---|
| 398 | *Afinal = A;
|
|---|
| 399 | *niter = 0;
|
|---|
| 400 | if(A<1.) return(-1.);
|
|---|
| 401 | if(A==1.) return(-2.);
|
|---|
| 402 | if(U<=0.) return(-3.);
|
|---|
| 403 |
|
|---|
| 404 | u1 = 0.;
|
|---|
| 405 | a1 = ampli_tf(u1,U,-1);
|
|---|
| 406 | if(a1==A) return(0.);
|
|---|
| 407 | /* on ne peut atteindre la valeur A */
|
|---|
| 408 | if(a1<A) { *Afinal = a1; return(-4.);}
|
|---|
| 409 |
|
|---|
| 410 | /* on recherche un intervalle encadrant A */
|
|---|
| 411 | u2 = u0_ponct(A);
|
|---|
| 412 | a2 = ampli_tf(u2,U,-1);
|
|---|
| 413 | while( a2>A && *niter<=nitermx ) {
|
|---|
| 414 | /* printf("recherche inter %f (%f), %f (%f)\n",u1,a1,u2,a2); */
|
|---|
| 415 | u1 = u2;
|
|---|
| 416 | a1 = a2;
|
|---|
| 417 | u2 *= 5.;
|
|---|
| 418 | a2 = ampli_tf(u2,U,-1);
|
|---|
| 419 | (*niter)++;
|
|---|
| 420 | }
|
|---|
| 421 |
|
|---|
| 422 | if( *niter>nitermx ) {
|
|---|
| 423 | *Afinal = a2;
|
|---|
| 424 | return(-5.);
|
|---|
| 425 | }
|
|---|
| 426 |
|
|---|
| 427 | /* a ce niveau u1 (a1) et u2 (a2) encadrent u0 (a0) */
|
|---|
| 428 | test = 100.*prec;
|
|---|
| 429 | *niter = 0;
|
|---|
| 430 | while (test > prec && *niter<= nitermx ) {
|
|---|
| 431 | /* printf("recherche conv %f (%f), %f (%f)\n",u1,a1,u2,a2); */
|
|---|
| 432 | u0 = (u1+u2)/2.;
|
|---|
| 433 | a0 = ampli_tf(u0,U,-1);
|
|---|
| 434 | (*niter)++;
|
|---|
| 435 | if(a0<A) {
|
|---|
| 436 | u2 = u0;
|
|---|
| 437 | a2 = a0;
|
|---|
| 438 | } else {
|
|---|
| 439 | u1 = u0;
|
|---|
| 440 | a1 = a0;
|
|---|
| 441 | }
|
|---|
| 442 | test = fabs(a0-A);
|
|---|
| 443 | }
|
|---|
| 444 |
|
|---|
| 445 | if( *niter > nitermx ) *niter = -1;
|
|---|
| 446 | /* printf("conv atteinte %f (%f), %f (%f)\n",u1,u2,a2); */
|
|---|
| 447 | *Afinal = a0;
|
|---|
| 448 | return(u0);
|
|---|
| 449 |
|
|---|
| 450 | }
|
|---|
| 451 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
