frechetslave.c

su 的源代码库

	    uuC2.r = 2 * uuC.r;
	    uuC2.i = 2 * uuC.i;
	    
	    muC.r = uC.r * -1;
	    muC.i = uC.i * -1;
	    
	    /* building reflectivity matrices */
	    RmFrechet();
	    Rp();
	 
	    /* reseting */
	    As1 = zeroC;      As2 = zeroC;      /* downgoing waves */
	    Cs1 = zeroC;      Cs2 = zeroC;      /* downgoing waves */
	    Bs1 = zeroC;      Bs2 = zeroC;      /* upgoing waves */
	    Ds1 = zeroC;      Ds2 = zeroC;      /* upgoing waves */
	    
	    /* P-wave potential */
	    /* PSlowness^2 - uuC */
	    auxm1 = PSlowness[0][0].r - uuC.r;
	    auxm2 = PSlowness[0][0].i - uuC.i;
	    auxm3 = sqrt(auxm1 * auxm1 + auxm2 * auxm2);
	    auxm3 = sqrt(auxm3);
	    angle = atan2(auxm2, auxm1) / 2;
	    am.r = auxm3 * cos(angle);
	    am.i = auxm3 * sin(angle);

	    /* am * I */
	    amI.r = -am.i;
	    amI.i = am.r;

	    As1 = uC;
	    if (directWave) Bs1 = muC;
	    
	    /* 1 / am */
	    aux = am.r * am.r + am.i * am.i;
	    amInv.r = am.r / aux;
	    amInv.i = -am.i / aux;
	    
	    /* amInv * uuC */
	    aux2.r = amInv.r * uuC.r - uuC.i * amInv.i;
	    aux2.i = amInv.r * uuC.i + amInv.i * uuC.r;

	    /* aux2 * -I */
	    As2.r = aux2.i;
	    As2.i = -aux2.r;

	    /* notice that Bs2 = As2 */
	    if (directWave) Bs2 = As2;
	    
	    /* S-wave potential */
	    /* SSlowness^2 - uuC */
	    auxm1 = SSlowness[0][0].r - uuC.r;
	    auxm2 = SSlowness[0][0].i - uuC.i;
	    
	    /* computing bm */
	    auxm3 = sqrt(auxm1 * auxm1 + auxm2 * auxm2);
	    auxm3 = sqrt(auxm3);
	    angle = atan2(auxm2, auxm1) / 2;
	    bm.r = auxm3 * cos(angle);
	    bm.i = auxm3 * sin(angle);
	    
	    /* bm * I */
	    bmI.r = -bm.i;
	    bmI.i = bm.r;

	    /* 1 / bm */
	    aux = bm.r * bm.r + bm.i * bm.i;
	    bmInv.r = bm.r / aux;
	    bmInv.i = -bm.i / aux;
	    
	    /* 1. / bm * uuC */
	    aux1.r = bmInv.r * uuC.r - bmInv.i * uuC.i;
	    aux1.i = bmInv.r * uuC.i + bmInv.i * uuC.r;
	 
	    /* notice that Cs1 = Ds1 */
	    Cs1 = aux1;
	    if (directWave) Ds1 = aux1;
	    
	    Cs2.r = -uC.i;
	    Cs2.i = uC.r;

	    if (directWave)
	    {
	       Ds2.r = -Cs2.r;
	       Ds2.i = -Cs2.i;
	    }

	    /* computing compensation for free-surface */
	    buildFreeSurfaceCompensation(am, bm);

	    /* computing phase shift (that's the matrix G in Muller's */
	    /* paper eq. (87) */
	    /* exp(j * am * wC * (-info->zs)) */
	    auxm1 = info->zs * (- amI.r * wC.r + amI.i * wC.i);
	    auxm2 = -info->zs * (amI.r * wC.i + amI.i * wC.r);
	    g[0].r = exp(auxm1) * cos(auxm2);
	    g[0].i = exp(auxm1) * sin(auxm2);

	    /* exp(j * bm * wC * (-info->zs)) */
	    auxm1 = info->zs * (- bmI.r * wC.r + bmI.i * wC.i);
	    auxm2 = -info->zs * (bmI.r * wC.i + bmI.i * wC.r);
	    g[1].r = exp(auxm1) * cos(auxm2);
	    g[1].i = exp(auxm1) * sin(auxm2);
	    
	    /* computing the product I - R-R+ */
	    auxm1 = rm[0][0].r * rp[0][0].r - rm[0][0].i * rp[0][0].i;
	    auxm2 = rm[0][0].r * rp[0][0].i + rm[0][0].i * rp[0][0].r;
	    auxm3 = rm[0][1].r * rp[1][0].r - rm[0][1].i * rp[1][0].i;
	    auxm4 = rm[0][1].r * rp[1][0].i + rm[0][1].i * rp[1][0].r;
	    irr[0][0].r = 1 - (auxm1 + auxm3);
	    irr[0][0].i = - (auxm2 + auxm4);

	    auxm1 = rm[0][0].r * rp[0][1].r - rm[0][0].i * rp[0][1].i;
	    auxm2 = rm[0][0].r * rp[0][1].i + rm[0][0].i * rp[0][1].r;
	    auxm3 = rm[0][1].r * rp[1][1].r - rm[0][1].i * rp[1][1].i;
	    auxm4 = rm[0][1].r * rp[1][1].i + rm[0][1].i * rp[1][1].r;
	    irr[0][1].r = - (auxm1 + auxm3);
	    irr[0][1].i = - (auxm2 + auxm4);
	    
	    auxm1 = rm[1][0].r * rp[0][0].r - rm[1][0].i * rp[0][0].i;
	    auxm2 = rm[1][0].r * rp[0][0].i + rm[1][0].i * rp[0][0].r;
	    auxm3 = rm[1][1].r * rp[1][0].r - rm[1][1].i * rp[1][0].i;
	    auxm4 = rm[1][1].r * rp[1][0].i + rm[1][1].i * rp[1][0].r;
	    irr[1][0].r = - (auxm1 + auxm3);
	    irr[1][0].i = - (auxm2 + auxm4);

	    auxm1 = rm[1][0].r * rp[0][1].r - rm[1][0].i * rp[0][1].i;
	    auxm2 = rm[1][0].r * rp[0][1].i + rm[1][0].i * rp[0][1].r;
	    auxm3 = rm[1][1].r * rp[1][1].r - rm[1][1].i * rp[1][1].i;
	    auxm4 = rm[1][1].r * rp[1][1].i + rm[1][1].i * rp[1][1].r;
	    irr[1][1].r = 1 - (auxm1 + auxm3);
	    irr[1][1].i = - (auxm2 + auxm4);

	    /* inverting irr explicitly */
	    auxm1 = irr[0][0].r * irr[1][1].r - irr[0][0].i * irr[1][1].i;
	    auxm2 = irr[0][0].r * irr[1][1].i + irr[0][0].i * irr[1][1].r;
	    auxm3 = irr[0][1].r * irr[1][0].r - irr[0][1].i * irr[1][0].i;
	    auxm4 = irr[0][1].r * irr[1][0].i + irr[0][1].i * irr[1][0].r;
	    aux1.r = auxm1 - auxm3;
	    aux1.i = auxm2 - auxm4;
	    
	    /* 1 / aux1 */
	    aux = aux1.r * aux1.r + aux1.i * aux1.i;
	    aux1.r = aux1.r / aux;
	    aux1.i = -aux1.i / aux;
	    
	    /* Inverse of irr */
	    irrI[0][0].r = irr[1][1].r * aux1.r - irr[1][1].i * aux1.i;
	    irrI[0][0].i = irr[1][1].r * aux1.i + irr[1][1].i * aux1.r;
	    
	    irrI[0][1].r = -(irr[0][1].r * aux1.r - irr[0][1].i * aux1.i);
	    irrI[0][1].i = -(irr[0][1].r * aux1.i + irr[0][1].i * aux1.r);
	    
	    irrI[1][0].r = -(irr[1][0].r * aux1.r - irr[1][0].i * aux1.i);
	    irrI[1][0].i = -(irr[1][0].r * aux1.i + irr[1][0].i * aux1.r);
	    
	    irrI[1][1].r = irr[0][0].r * aux1.r - irr[0][0].i * aux1.i;
	    irrI[1][1].i = irr[0][0].r * aux1.i + irr[0][0].i * aux1.r;
	    
	    /* computing vectors V1,2, check eq (76) Muller's paper */
	    auxm1 = As1.r * rm[0][0].r - As1.i * rm[0][0].i;
	    auxm2 = As1.r * rm[0][0].i + As1.i * rm[0][0].r;
	    auxm3 = Cs1.r * rm[0][1].r - Cs1.i * rm[0][1].i;
	    auxm4 = Cs1.r * rm[0][1].i + Cs1.i * rm[0][1].r;
	    aux1.r = Bs1.r + (auxm1 + auxm3);
	    aux1.i = Bs1.i + (auxm2 + auxm4);
	    
	    auxm1 = As1.r * rm[1][0].r - As1.i * rm[1][0].i;
	    auxm2 = As1.r * rm[1][0].i + As1.i * rm[1][0].r;
	    auxm3 = Cs1.r * rm[1][1].r - Cs1.i * rm[1][1].i;
	    auxm4 = Cs1.r * rm[1][1].i + Cs1.i * rm[1][1].r;
	    aux2.r = Ds1.r + (auxm1 + auxm3);
	    aux2.i = Ds1.i + (auxm2 + auxm4);
	    
	    auxm1 = aux1.r * irrI[0][0].r - aux1.i * irrI[0][0].i;
	    auxm2 = aux1.r * irrI[0][0].i + aux1.i * irrI[0][0].r;
	    auxm3 = aux2.r * irrI[0][1].r - aux2.i * irrI[0][1].i;
	    auxm4 = aux2.r * irrI[0][1].i + aux2.i * irrI[0][1].r;
	    v1[0][0].r = auxm1 + auxm3;
	    v1[0][0].i = auxm2 + auxm4;
	    
	    auxm1 = aux1.r * irrI[1][0].r - aux1.i * irrI[1][0].i;
	    auxm2 = aux1.r * irrI[1][0].i + aux1.i * irrI[1][0].r;
	    auxm3 = aux2.r * irrI[1][1].r - aux2.i * irrI[1][1].i;
	    auxm4 = aux2.r * irrI[1][1].i + aux2.i * irrI[1][1].r;
	    v1[0][1].r = auxm1 + auxm3;
	    v1[0][1].i = auxm2 + auxm4;
	    
	    /* loop over "active" layers */
	    for (iDer = 1, i = 0; i < numberPar; i++)
	    {
	       /* i = 0 -> Vp  */
	       /* i = 1 -> Vs  */
	       /* i = 2 -> rho */
	       for (iL = MIN(lim[0], 2); iL < MIN(lim[0], 2) + limRange; 
		    iL++, iDer++)
	       {
		  /* rp * [v1[0], v1[1]] + (As1, Cs1)*/
		  auxm1 = rp[0][0].r * v1[0][0].r - rp[0][0].i * v1[0][0].i;
		  auxm2 = rp[0][0].r * v1[0][0].i + rp[0][0].i * v1[0][0].r;
		  auxm1 += rp[0][1].r * v1[0][1].r - rp[0][1].i * v1[0][1].i 
	  	        + As1.r;
		  auxm2 += rp[0][1].r * v1[0][1].i + rp[0][1].i * v1[0][1].r 
		        + As1.i;
	    
		  auxm3 = rp[1][0].r * v1[0][0].r - rp[1][0].i * v1[0][0].i;
		  auxm4 = rp[1][0].r * v1[0][0].i + rp[1][0].i * v1[0][0].r;
		  auxm3 += rp[1][1].r * v1[0][1].r - rp[1][1].i * v1[0][1].i 
		        + Cs1.r;
		  auxm4 += rp[1][1].r * v1[0][1].i + rp[1][1].i * v1[0][1].r 
		        + Cs1.i;

		  /* DmB[0][active layers][0 1 2 3] * */
		  /*              ((auxm1, auxm2), (auxm3, auxm4)) */
		  aux1.r = auxm1 * DmB[0][i * limRange + iL][0].r 
		         - auxm2 * DmB[0][i * limRange + iL][0].i 
			 + auxm3 * DmB[0][i * limRange + iL][1].r 
			 - auxm4 * DmB[0][i * limRange + iL][1].i;
		  aux1.i = auxm1 * DmB[0][i * limRange + iL][0].i 
		         + auxm2 * DmB[0][i * limRange + iL][0].r
			 + auxm3 * DmB[0][i * limRange + iL][1].i 
			 + auxm4 * DmB[0][i * limRange + iL][1].r;
	       
		  aux2.r = auxm1 * DmB[0][i * limRange + iL][2].r 
		         - auxm2 * DmB[0][i * limRange + iL][2].i 
		         + auxm3 * DmB[0][i * limRange + iL][3].r 
		         - auxm4 * DmB[0][i * limRange + iL][3].i;
		  aux2.i = auxm1 * DmB[0][i * limRange + iL][2].i 
		         + auxm2 * DmB[0][i * limRange + iL][2].r  
		         + auxm3 * DmB[0][i * limRange + iL][3].i 
		         + auxm4 * DmB[0][i * limRange + iL][3].r;
	       
		  /* irrI * (aux1, aux2) */
		  auxm1 = irrI[0][0].r * aux1.r - irrI[0][0].i * aux1.i;
		  auxm2 = irrI[0][0].r * aux1.i + irrI[0][0].i * aux1.r;
		  auxm3 = irrI[0][1].r * aux2.r - irrI[0][1].i * aux2.i;
		  auxm4 = irrI[0][1].r * aux2.i + irrI[0][1].i * aux2.r;
		  v1[iDer][0].r = auxm1 + auxm3;
		  v1[iDer][0].i = auxm2 + auxm4;
		  
		  auxm1 = irrI[1][0].r * aux1.r - irrI[1][0].i * aux1.i;
		  auxm2 = irrI[1][0].r * aux1.i + irrI[1][0].i * aux1.r;
		  auxm3 = irrI[1][1].r * aux2.r - irrI[1][1].i * aux2.i;
		  auxm4 = irrI[1][1].r * aux2.i + irrI[1][1].i * aux2.r;
		  v1[iDer][1].r = auxm1 + auxm3;
		  v1[iDer][1].i = auxm2 + auxm4;
		  
		  /* DD 
		     fprintf(stderr, 
		     "iDer %d i %d iL % d DMB[[%f %f %f %f]\n", 
		     iDer, i, iL, DmB[0][i * limRange + iL][0].r, 
		     DmB[0][i * limRange + iL][1].r, 
		     DmB[0][i * limRange + iL][2].r,
		     DmB[0][i * limRange + iL][3].r);*/
	       }
	    }
	    
	    auxm1 = As2.r * rm[0][0].r - As2.i * rm[0][0].i;
	    auxm2 = As2.r * rm[0][0].i + As2.i * rm[0][0].r;
	    auxm3 = Cs2.r * rm[0][1].r - Cs2.i * rm[0][1].i;
	    auxm4 = Cs2.r * rm[0][1].i + Cs2.i * rm[0][1].r;
	    aux1.r = Bs2.r + (auxm1 + auxm3);
	    aux1.i = Bs2.i + (auxm2 + auxm4);
	    
	    auxm1 = As2.r * rm[1][0].r - As2.i * rm[1][0].i;
	    auxm2 = As2.r * rm[1][0].i + As2.i * rm[1][0].r;
	    auxm3 = Cs2.r * rm[1][1].r - Cs2.i * rm[1][1].i;
	    auxm4 = Cs2.r * rm[1][1].i + Cs2.i * rm[1][1].r;
	    aux2.r = Ds2.r + (auxm1 + auxm3);
	    aux2.i = Ds2.i + (auxm2 + auxm4);
	    
	    auxm1 = aux1.r * irrI[0][0].r - aux1.i * irrI[0][0].i;
	    auxm2 = aux1.r * irrI[0][0].i + aux1.i * irrI[0][0].r;
	    auxm3 = aux2.r * irrI[0][1].r - aux2.i * irrI[0][1].i;
	    auxm4 = aux2.r * irrI[0][1].i + aux2.i * irrI[0][1].r;
	    v2[0][0].r = auxm1 + auxm3;
	    v2[0][0].i = auxm2 + auxm4;
	    
	    auxm1 = aux1.r * irrI[1][0].r - aux1.i * irrI[1][0].i;
	    auxm2 = aux1.r * irrI[1][0].i + aux1.i * irrI[1][0].r;
	    auxm3 = aux2.r * irrI[1][1].r - aux2.i * irrI[1][1].i;
	    auxm4 = aux2.r * irrI[1][1].i + aux2.i * irrI[1][1].r;
	    v2[0][1].r = auxm1 + auxm3;
	    v2[0][1].i = auxm2 + auxm4;

	    /* loop over "active" layers */
	    for (iDer = 1, i = 0; i < numberPar; i++)
	    {
	       /* i = 0 -> Vp  */
	       /* i = 1 -> Vs  */
	       /* i = 2 -> rho */
	       for (iL = MIN(lim[0], 2); iL < MIN(lim[0], 2) + limRange; 
		    iL++, iDer++)
	       {
		  /* rp * [v2[0], v2[1]] + (As2, Bs2) */
		  auxm1 = rp[0][0].r * v2[0][0].r - rp[0][0].i * v2[0][0].i;
		  auxm2 = rp[0][0].r * v2[0][0].i + rp[0][0].i * v2[0][0].r;
		  auxm1 += rp[0][1].r * v2[0][1].r - rp[0][1].i * v2[0][1].i 
	 	        + As2.r;
		  auxm2 += rp[0][1].r * v2[0][1].i + rp[0][1].i * v2[0][1].r 
		        + As2.i;

		  auxm3 = rp[1][0].r * v2[0][0].r - rp[1][0].i * v2[0][0].i;
		  auxm4 = rp[1][0].r * v2[0][0].i + rp[1][0].i * v2[0][0].r;
		  auxm3 += rp[1][1].r * v2[0][1].r - rp[1][1].i * v2[0][1].i 
	 	        + Cs2.r;
		  auxm4 += rp[1][1].r * v2[0][1].i + rp[1][1].i * v2[0][1].r 
		        + Cs2.i;

		  /* DmB[0][active layers][0 1 2 3] * */
		  /*                      ((auxm1, auxm2), (auxm3, auxm4)) */
		  aux1.r = auxm1 * DmB[0][i * limRange + iL][0].r 
		         - auxm2 * DmB[0][i * limRange + iL][0].i 
			 + auxm3 * DmB[0][i * limRange + iL][1].r 
			 - auxm4 * DmB[0][i * limRange + iL][1].i;
		  aux1.i = auxm1 * DmB[0][i * limRange + iL][0].i 
		         + auxm2 * DmB[0][i * limRange + iL][0].r 
			 + auxm3 * DmB[0][i * limRange + iL][1].i 
			 + auxm4 * DmB[0][i * limRange + iL][1].r;
		  aux2.r = auxm1 * DmB[0][i * limRange + iL][2].r 
		         - auxm2 * DmB[0][i * limRange + iL][2].i  
			 + auxm3 * DmB[0][i * limRange + iL][3].r 
			 - auxm4 * DmB[0][i * limRange + iL][3].i;
		  aux2.i = auxm1 * DmB[0][i * limRange + iL][2].i 
		         + auxm2 * DmB[0][i * limRange + iL][2].r 
			 + auxm3 * DmB[0][i * limRange + iL][3].i 
			 + auxm4 * DmB[0][i * limRange + iL][3].r;