shreflect.c

su 的源代码库

							td0n22[ip]=cmplx(0.0,0.0);
							ewtu211=cmul(rvrb211,cmul(ewigh11,tu11));
							tun0p11=cmul(ewtu211,tun011[ip]);
							tun0p21=cmul(ewtu211,tun021[ip]);
							rtb111=cmul(ewrd11,ewd0211);
							rtb112=cmul(ewrd11,ewd0212);
							rd0np11=cadd(rd0n11[ip],cmul(tun011[ip],rtb111));
							rd0np12=cadd(rd0n12[ip],cmul(tun011[ip],rtb112));
							rd0np21=cadd(rd0n21[ip],cmul(tun021[ip],rtb111));
							rd0np22=cadd(rd0n22[ip],cmul(tun022[ip],rtb112));
							run011[ip]=cadd(ru11,cmul(ewtu211,
								cmul(wrn011,td11)));
							run012[ip]=cmplx(0.0,0.0);
							run021[ip]=cmplx(0.0,0.0);
							run022[ip]=cmplx(0.0,0.0);
							tun011[ip]=tun0p11;
							tun012[ip]=cmplx(0.0,0.0);
							tun021[ip]=tun0p21;
							tun022[ip]=cmplx(0.0,0.0);
							rd0n11[ip]=rd0np11;
							rd0n12[ip]=rd0np12;
							rd0n21[ip]=rd0np21;
							rd0n22[ip]=rd0np22;
						}
					}
				} 
				/* if (il==lobs[ijk]-1) { */
				if (il==lobs[ijk]-2) {
					ik1=ijk*block_size;		
					for (ip=0; ip<block_size; ip++) {
						iz=ik1+ip;
						rurf11[iz]=run011[ip];
						rurf12[iz]=run012[ip];
						rurf21[iz]=run021[ip];
						rurf22[iz]=run022[ip];
						rdfr11[iz]=rd0n11[ip];
						rdfr12[iz]=rd0n12[ip];
						rdfr21[iz]=rd0n21[ip];
						rdfr22[iz]=rd0n22[ip];
						turf11[iz]=tun011[ip];
						turf12[iz]=tun012[ip];
						turf21[iz]=tun021[ip];
						turf22[iz]=tun022[ip];
						tdfr11[iz]=td0n11[ip];
						tdfr12[iz]=td0n12[ip];
						tdfr21[iz]=td0n21[ip];
						tdfr22[iz]=td0n22[ip];
					}
					ijk++;
				} 

				/* if (il==lsource-1) { */
				if (il==lsource-2) {
			    for (ip=0; ip<block_size; ip++) {
				rusf11[ip]=run011[ip];
				rusf12[ip]=run012[ip];
				rusf21[ip]=run021[ip];
				rusf22[ip]=run022[ip];
				rdfs11[ip]=rd0n11[ip];
				rdfs12[ip]=rd0n12[ip];
				rdfs21[ip]=rd0n21[ip];
				rdfs22[ip]=rd0n22[ip];
				run011[ip]=cmplx(0.0,0.0);
				run012[ip]=cmplx(0.0,0.0);
				run021[ip]=cmplx(0.0,0.0);
				run022[ip]=cmplx(0.0,0.0);
				rd0n11[ip]=cmplx(0.0,0.0);
				rd0n12[ip]=cmplx(0.0,0.0);
				rd0n21[ip]=cmplx(0.0,0.0);
				rd0n22[ip]=cmplx(0.0,0.0);
				tusf11[ip]=tun011[ip];
				tusf12[ip]=tun012[ip];
				tusf21[ip]=tun021[ip];
				tusf22[ip]=tun022[ip];
				tun011[ip]=cmplx(1.0,0.0);
				tun012[ip]=cmplx(0.0,0.0);
				tun021[ip]=cmplx(0.0,0.0);
				tun022[ip]=cmplx(1.0,0.0);
				td0n11[ip]=cmplx(1.0,0.0);
				td0n12[ip]=cmplx(0.0,0.0);
				td0n21[ip]=cmplx(0.0,0.0);
				   	td0n22[ip]=cmplx(1.0,0.0);
	    		}
			}
    		}

			/* loop over receiver depths */
		    for (iz=0; iz<nor; iz++) {
			ik1=iz*block_size;
			for (ip=0; ip<block_size; ip++) {
			    ijk1=ik1+ip;    
			    det=csub(rd0n11[ip],rdfr11[ijk1]);
			    if (rcabs(det)<=SZERO) {
				prv[iz]=flag[iz];
				flag[iz]=2;
				ibl[iz]=ip;
				break;
			    }
			}
    		}


 		   for (iz=0; iz<nor; iz++) {
		ik1=iz*block_size;
			lrec=lobs[iz];
			ik2=(lrec-1)*block_size;

			/* receiver above the source */
			if (lsource>lrec) {
			    for (ip=0; ip<block_size; ip++) {
				ijk1=ik1+ip;
						ijk2=ik2+ip;
						rtb111=cdiv(cmplx(1.0,0.0),turf11[ijk1]);
						rdnr11[ip]=cmul(rtb111,tusf11[ip]);
						rdnr12[ip]=cmul(rtb111,tusf12[ip]);
						rdnr21[ip]=cmplx(0.0,0.0);
						rdnr22[ip]=cmplx(0.0,0.0);
						rtb111=rdnr11[ip];
						rtb112=rdnr12[ip];
						rtb121=rdnr21[ip];
						rtb122=rdnr22[ip];
						rtb211=cadd(cmul(rurf11[ijk1],rtb111),cmul(rurf12[ijk1],
							rtb121));
						rtb212=cadd(cmul(rurf11[ijk1],rtb112),cmul(rurf12[ijk1],
							rtb122));
						rtb221=cadd(cmul(rurf21[ijk1],rtb111),cmul(rurf22[ijk1],
							rtb121));
						rtb222=cadd(cmul(rurf21[ijk1],rtb112),cmul(rurf22[ijk1],
							rtb122));
						rtb311=cadd(cmul(rd0n11[ip],rusf11[ip]),cmul(rd0n12[ip],
							rusf21[ip]));
						rtb312=cadd(cmul(rd0n11[ip],rusf12[ip]),cmul(rd0n12[ip],
							rusf22[ip]));
						rtb321=cadd(cmul(rd0n21[ip],rusf11[ip]),cmul(rd0n22[ip],
							rusf21[ip]));
						rtb322=cadd(cmul(rd0n21[ip],rusf12[ip]),cmul(rd0n22[ip],
							rusf22[ip]));
						x=csub(cmplx(1.0,0.0),rtb311);
						y=csub(cmplx(1.0,0.0),rtb322);	
						det=csub(cmul(x,y),cmul(rtb312,rtb321));
						rvrb111=cdiv(y,det);
						rvrb112=cdiv(rtb312,det);
						rvrb121=cdiv(rtb321,det);
						rvrb122=cdiv(x,det);
						fact11=cadd(cmul(rd0n11[ip],sigmad1[ip]),
							cadd(cmul(rd0n12[ip],sigmad2[ip]),sigmau1[ip]));
						fact12=cadd(cmul(rd0n21[ip],sigmad1[ip]),
							cadd(cmul(rd0n22[ip],sigmad2[ip]),sigmau2[ip]));
						vuzs1=cadd(cmul(rvrb111,fact11),cmul(rvrb112,fact12));
						vuzs2=cadd(cmul(rvrb121,fact11),cmul(rvrb122,fact12));
						up[ip]=cmul(p[ip],cadd(cmul(cadd(rtb111,rtb211),vuzs1),
							cmul(cadd(rtb112,rtb212),vuzs2)));
					}
				} else {
				
					/* receiver below the source */
					if (flag[iz]==2) {
						for (ip=0; ip<ibl[iz]-1; ip++) {
							ijk1=ik1+ip;
							rtb211=cdiv(turf11[ijk1],csub(rd0n11[ip],
							rdfr11[ijk1]));
							rdnr11[ip]=cdiv(cmplx(1.0,0.0),cadd(cmul(tdfr11[ijk1],rtb211),
								rurf11[ijk1]));
							rdnr12[ip]=cmplx(0.0,0.0);   
							rdnr21[ip]=cmplx(0.0,0.0);   
							rdnr22[ip]=cmplx(0.0,0.0);   
						}
						for (ip=ibl[iz]-1; ip<block_size; ip++) {
							rdnr11[ip]=cmplx(0.0,0.0);
							rdnr12[ip]=cmplx(0.0,0.0);
							rdnr21[ip]=cmplx(0.0,0.0);
							rdnr22[ip]=cmplx(0.0,0.0);
						}
					} else if (flag[iz]==0) {
			
						for (ip=0; ip<block_size; ip++) {
							ijk1=ik1+ip;
							rtb211=cdiv(turf11[ijk1],csub(rd0n11[ip],
								rdfr11[ijk1]));
							rdnr11[ip]=cdiv(cmplx(1.0,0.0),cadd(cmul(tdfr11[ijk1],rtb211),
								rurf11[ijk1]));
							rdnr12[ip]=cmplx(0.0,0.0);
							rdnr21[ip]=cmplx(0.0,0.0);
							rdnr22[ip]=cmplx(0.0,0.0);
						}
					} else {
						for (ip=0; ip<block_size; ip++) {
							rdnr11[ip]=cmplx(0.0,0.0);
							rdnr12[ip]=cmplx(0.0,0.0);
							rdnr21[ip]=cmplx(0.0,0.0);
							rdnr22[ip]=cmplx(0.0,0.0);
						}
					}	
					if (flag[iz]==2) flag[iz]=prv[iz];	
					for (ip=0; ip<block_size; ip++) {
						ijk1=ik1+ip;
						ijk2=ik2+ip;
						rtb111=cadd(cmul(rurf11[ijk1],rdnr11[ip]),
							cmul(rurf12[ijk1],rdnr21[ip]));
						rtb112=cadd(cmul(rurf11[ijk1],rdnr12[ip]),
							cmul(rurf12[ijk1],rdnr22[ip]));
						rtb121=cadd(cmul(rurf21[ijk1],rdnr11[ip]),
							cmul(rurf22[ijk1],rdnr21[ip]));
						rtb122=cadd(cmul(rurf21[ijk1],rdnr12[ip]),
							cmul(rurf22[ijk1],rdnr22[ip]));
						x=csub(cmplx(1.0,0.0),rtb111);
						y=csub(cmplx(1.0,0.0),rtb122);
						det=csub(cmul(x,y),cmul(rtb112,rtb121));
						rvrb111=cdiv(y,det);
						rvrb112=cdiv(rtb112,det);
						rvrb121=cdiv(rtb121,det);
						rvrb122=cdiv(x,det);
						rtb111=cadd(cmul(rvrb111,tdfr11[ijk1]),cmul(rvrb112,
							tdfr21[ijk1]));
						rtb112=cadd(cmul(rvrb111,tdfr12[ijk1]),cmul(rvrb112,
							tdfr22[ijk1])); 
						rtb121=cadd(cmul(rvrb121,tdfr11[ijk1]),cmul(rvrb122,
							tdfr21[ijk1]));
						rtb122=cadd(cmul(rvrb121,tdfr12[ijk1]),cmul(rvrb122,
							tdfr22[ijk1]));
						rtb211=cadd(cmul(rdnr11[ip],rtb111),cmul(rdnr12[ip],
							rtb121));
						rtb212=cadd(cmul(rdnr11[ip],rtb112),cmul(rdnr12[ip],
							rtb122));
						rtb221=cadd(cmul(rdnr21[ip],rtb111),cmul(rdnr22[ip],
							rtb121));
						rtb222=cadd(cmul(rdnr21[ip],rtb112),cmul(rdnr22[ip],
							rtb122));
						rtb311=cadd(cmul(rusf11[ip],rd0n11[ip]),cmul(rusf12[ip],
							rd0n21[ip]));
						rtb312=cadd(cmul(rusf11[ip],rd0n12[ip]),cmul(rusf12[ip],
							rd0n22[ip]));
						rtb321=cadd(cmul(rusf21[ip],rd0n11[ip]),cmul(rusf22[ip],
							rd0n21[ip]));
						rtb322=cadd(cmul(rusf21[ip],rd0n12[ip]),cmul(rusf22[ip],
							rd0n22[ip]));
						x=csub(cmplx(1.0,0.0),rtb311);
						y=csub(cmplx(1.0,0.0),rtb322);
						det=csub(cmul(x,y),cmul(rtb312,rtb321));
						rvrb111=cdiv(y,det);
						rvrb112=cdiv(rtb312,det);
						rvrb121=cdiv(rtb321,det);
						rvrb122=cdiv(x,det);
						fact11=cadd(cmul(rusf11[ip],sigmau1[ip]),cadd(
							cmul(rusf12[ip],sigmau2[ip]),sigmad1[ip]));
						fact12=cadd(cmul(rusf21[ip],sigmau1[ip]),cadd(
							cmul(rusf22[ip],sigmau2[ip]),sigmad2[ip]));
						vdzs1=cadd(cmul(rvrb111,fact11),cmul(rvrb112,fact12));
						vdzs2=cadd(cmul(rvrb121,fact11),cmul(rvrb122,fact12));
						up[ip]=cmul(p[ip],cadd(cmul(cadd(rtb111,rtb211),vdzs1),
							cmul(cadd(rtb112,rtb212),vdzs2)));
					}
				}
				for (ix=0; ix<nx; ix++) {
					x1=bx+dx*ix;
					if (x1==0.0) {
						for (ip=1; ip<block_size-1; ip++) {
							t1=cdiv(cmplx(1.0,0.0),csqrt(p[ip]));
							vos1[ip]=cmul(up[ip],t1);
						}
						vos1[0]=cmplx(0.0,0.0);
						vos1[1]=crmul(vos1[1],0.5);
						vos1[block_size-2]=crmul(vos1[block_size-2],.5);
					} else {
						temr=crmul(wpie,2*PI*x1);
						temr=cdiv(tem,csqrt(temr));
						sig=crmul(divfac,x1);
						sigh=cmul(sig,cdp);
						for (ip=0; ip<block_size-1; ip++) {
							texp[ip]=cexp(cmul(sig,p[ip]));
						}

						/******************************************************
						*			Compute the slowness integral			  *
						******************************************************/
						if (int_type==1) {		/* use trapezoidal rule */
							temr=crmul(cmul(temr,cdp),0.5);
							for (ip=0; ip<block_size-1; ip++) {
								jj=ip+1;
								vos1[ip]=cmul(temr,cadd(cmul(up[ip],texp[ip]),
									cmul(up[jj],texp[jj])));
							}
						} else if (int_type==2) {
							for (ip=0; ip<block_size-1; ip++) {
								jj=ip+1;
								t1=csub(texp[jj],texp[ip]);
								func1=csub(cmul(up[jj],texp[jj]),
									cmul(up[ip],texp[ip]));
								func2=cmul(csub(up[jj],up[ip]),t1);
								vos1[ip]=cmul(cdiv(cdp,sigh),cmul(csub(func1,
									cdiv(func2,sigh)),temr));
							}
						}
					}
		/**********************************************************
	    *		update output array			*
	    **********************************************************/
				for (ip=0; ip<block_size-1; ip++) {
						response1[iw][ix][iz]=cadd(response1[iw][ix][iz],
							vos1[ip]);
					}
				}
			}

	/* update loop variables */
            nblock--;
	    if (nblock != 0) {

	/*	bp=p[block_size].r; */
		bp=p[block_size-1].r;
		left++;
		if (left > block_size) {
		    left -=block_size;
		    nblock++;
		}

	    } else {
				if (int_type==1) bp=bp1;
				else if (int_type==2) bp=0.0;
				break;	/* last block, exit while loop */
			}
		}
	}
	
	/* free allocated space */
	free1complex(up);
	free1complex(rd0n11);free1complex(rd0n12);free1complex(rd0n21);
	free1complex(rd0n22);free1complex(td0n11);free1complex(td0n12);
	free1complex(td0n22);free1complex(tun011);free1complex(tun012);
	free1complex(tun021);free1complex(tun022);free1complex(run011);
	free1complex(run012);free1complex(run021);free1complex(run022);
	free1complex(tusf11);free1complex(tusf12);free1complex(tusf21);
	free1complex(tusf22);free1complex(rdfs11);free1complex(rdfs12);
	free1complex(rdfs21);free1complex(rdfs22);free1complex(rurf11);
	free1complex(rurf12);free1complex(rurf21);free1complex(rurf22);
	free1complex(turf11);free1complex(turf12);free1complex(turf21);
	free1complex(turf22);free1complex(rdfr11);free1complex(rdfr12);
	free1complex(rdfr21);free1complex(rdfr22);free1complex(tdfr11);
	free1complex(tdfr12);free1complex(tdfr21);free1complex(tdfr22);
	free1complex(r11);free1complex(r12);free1complex(r21);free1complex(r22);
	free1complex(vos1);free1complex(vos2);free1complex(vos3);

	/* free working space */
    free1complex(p);
    free1complex(pp);
    free1complex(pwin);
    free1complex(gl);
    free1complex(gt);
    free1complex(gam);
    free1complex(alpha);
    free1complex(betha);
    free1complex(sigmau1);
    free1complex(sigmau2);
    free1complex(sigmad1);
    free1complex(sigmad2);

}