psvreflect.c

su 的源代码库

							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);
							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);
						}
					} else {
						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[ip]=rd0n11[ip];
							rdfr12[ip]=rd0n12[ip];
							rdfr21[ip]=rd0n21[ip];
							rdfr22[ip]=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-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];
						tdfs11[ip]=td0n11[ip];
						tdfs12[ip]=td0n12[ip];
						tdfs21[ip]=td0n21[ip];
						tdfs22[ip]=td0n22[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 receivers */
			flg1=0;
			for (iz=0; iz<nor; iz++) {
				ik1=iz*block_size;
				for (ip=0; ip<block_size; ip++) {
					ijk1=ik1+ip;	
					if (acoustic[iz]==1) {
						det=csub(rd0n11[ip],rdfr11[ijk1]);
					} else {
						det=csub(cmul(csub(rd0n11[ip],rdfr11[ijk1]),
							csub(rd0n22[ip],rdfr22[ijk1])),cmul(csub(rd0n12[ip],
							rdfr12[ijk1]),csub(rd0n21[ip],rdfr21[ijk1])));
					}	
					if (rcabs(det)<=SZERO) {
						prv[iz]=flag[iz];
						flag[iz]=2;
						ibl[iz]=ip;
						break;
					}
				}
				if (acoustic[iz]==2) {
					flg1=1;
					for (ip=0; ip<block_size; ip++) {
						ijk1=ik1+ip;
						rnr11[ip]=rd0n11[ip];
						rnr12[ip]=rd0n12[ip];
						rnr21[ip]=rd0n21[ip];
						rnr22[ip]=rd0n22[ip];
						tdsr11[ip]=tdfr11[ijk1];
						tdsr12[ip]=tdfr12[ijk1];
						tdsr21[ip]=tdfr21[ijk1];
						tdsr22[ip]=tdfr22[ijk1];
						tusr11[ip]=turf11[ijk1];
						tusr12[ip]=turf12[ijk1];
						tusr21[ip]=turf21[ijk1];
						tusr22[ip]=turf22[ijk1];
						rdsr11[ip]=rdfr11[ijk1];	
						rdsr12[ip]=rdfr12[ijk1];	
						rdsr21[ip]=rdfr21[ijk1];	
						rdsr22[ip]=rdfr22[ijk1];	
						rusr11[ip]=rurf11[ijk1];
						rusr12[ip]=rurf12[ijk1];
						rusr21[ip]=rurf21[ijk1];
						rusr22[ip]=rurf22[ijk1];
						rtb111=cadd(cmul(rurf11[ijk1],rnr11[ip]),
							cmul(rurf12[ijk1],rnr21[ip]));
						rtb112=cadd(cmul(rurf11[ijk1],rnr12[ip]),
							cmul(rurf12[ijk1],rnr22[ip]));
						rtb121=cadd(cmul(rurf21[ijk1],rnr11[ip]),
							cmul(rurf22[ijk1],rnr21[ip]));
						rtb122=cadd(cmul(rurf21[ijk1],rnr12[ip]),
							cmul(rurf22[ijk1],rnr22[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(rnr11[ip],rtb111),
							cmul(rnr12[ip],rtb121));	
						rtb212=cadd(cmul(rnr11[ip],rtb112),
							cmul(rnr12[ip],rtb122));	
						rtb221=cadd(cmul(rnr21[ip],rtb111),
							cmul(rnr22[ip],rtb121));	
						rtb222=cadd(cmul(rnr21[ip],rtb112),
							cmul(rnr22[ip],rtb122));	
						rd0n11[ip]=cadd(cmul(rtb211,turf11[ijk1]),
							cadd(cmul(rtb221,turf12[ijk1]),rdfr11[ijk1]));
						rd0n12[ip]=cadd(cmul(rtb212,turf11[ijk1]),
							cadd(cmul(rtb222,turf12[ijk1]),rdfr12[ijk1]));
						rd0n21[ip]=cadd(cmul(rtb211,turf21[ijk1]),
							cadd(cmul(rtb221,turf22[ijk1]),rdfr21[ijk1]));
						rd0n22[ip]=cadd(cmul(rtb212,turf21[ijk1]),
							cadd(cmul(rtb222,turf22[ijk1]),rdfr22[ijk1]));
					}
				}
			}

			/* loop over receivers */
			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) {
					if (acoustic[iz]==4) {
						for (ip=0; ip<block_size; ip++) {
							ijk1=ik1+ip;
							ijk2=ik2+ip;
							rtb111=cadd(cmul(rdfs11[ip],rurf11[ijk1]),
								cmul(rdfs12[ip],rurf21[ijk1]));
							rtb112=cadd(cmul(rdfs11[ip],rurf12[ijk1]),
								cmul(rdfs12[ip],rurf22[ijk1]));
							rtb121=cadd(cmul(rdfs21[ip],rurf11[ijk1]),
								cmul(rdfs22[ip],rurf21[ijk1]));
							rtb122=cadd(cmul(rdfs21[ip],rurf12[ijk1]),
								cmul(rdfs22[ip],rurf22[ijk1]));
							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,tusf11[ip]),cmul(rvrb112,
								tusf21[ip]));
							rtb112=cadd(cmul(rvrb111,tusf12[ip]),cmul(rvrb112,
								tusf22[ip]));
							rtb121=cadd(cmul(rvrb121,tusf11[ip]),cmul(rvrb122,
								tusf21[ip]));
							rtb122=cadd(cmul(rvrb121,tusf12[ip]),cmul(rvrb122,
								tusf22[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(tdfs11[ip],rtb211),cmul(tdfs12[ip],
								rtb221));
							rtb312=cadd(cmul(tdfs11[ip],rtb212),cmul(tdfs12[ip],
								rtb222));
							rtb321=cadd(cmul(tdfs21[ip],rtb211),cmul(tdfs22[ip],
								rtb221));
							rtb322=cadd(cmul(tdfs21[ip],rtb212),cmul(tdfs22[ip],
								rtb222));
							rusf11[ip]=cadd(rusf11[ip],rtb311);
							rusf12[ip]=cadd(rusf12[ip],rtb312);
							rusf21[ip]=cadd(rusf21[ip],rtb321);
							rusf22[ip]=cadd(rusf22[ip],rtb322);
							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(sigmad1[ip],rd0n11[ip]),
								cadd(sigmau1[ip],cmul(rd0n12[ip],sigmad2[ip])));
							fact12=cadd(cmul(rd0n21[ip],sigmad1[ip]),
								cadd(sigmau2[ip],cmul(rd0n22[ip],sigmad2[ip])));
							vuzs1=cadd(cmul(rvrb111,fact11),
								cmul(rvrb112,fact12));
							vuzs2=cadd(cmul(rvrb121,fact11),
								cmul(rvrb122,fact12));
							up[ip]=cmul(prs[iz],cadd(cmul(vuzs1,cadd(rtb111,
								rtb211)),cmul(vuzs2,cadd(rtb112,rtb212))));
							rtb311=cadd(cmul(rurf11[ijk1],p[ip]),
								csub(cmul(p[ip],gt[ijk2]),rurf21[ijk1]));
							rtb312=cadd(cmul(rurf12[ijk1],p[ip]),
								csub(cmul(p[ip],gt[ijk2]),rurf22[ijk1]));
							rtb321=cadd(cmul(rurf11[ijk1],gl[ijk2]),
								csub(cmul(p[ip],rurf21[ijk1]),gl[ijk2]));
							rtb322=cadd(cmul(rurf12[ijk1],gl[ijk2]),
								cadd(cmul(p[ip],rurf22[ijk1]),p[ip]));
							rtb211=cadd(cmul(rtb311,rtb111),
								cmul(rtb312,rtb121));
							rtb212=cadd(cmul(rtb311,rtb112),
								cmul(rtb312,rtb122));
							rtb221=cadd(cmul(rtb321,rtb111),
								cmul(rtb322,rtb121));
							rtb222=cadd(cmul(rtb321,rtb112),
								cmul(rtb322,rtb122));
							ux[ip]=cadd(cmul(rtb211,vuzs1),cmul(rtb212,vuzs2));
							uz[ip]=cadd(cmul(rtb221,vuzs1),cmul(rtb222,vuzs2));
						}
					} else {
	
						/* source elastic,receiver elastic or*/
						/* source acoustic, receiver elastic */
						if (acoustic[iz]==0||acoustic[iz]==3) {
							for (ip=0; ip<block_size; ip++) {
								ijk1=ik1+ip;
								det=csub(cmul(turf11[ijk1],turf22[ijk1]),
									cmul(turf21[ijk1],turf12[ijk1]));
								rtb211=cdiv(turf22[ijk1],det);
								rtb212=cneg(cdiv(turf12[ijk1],det));
								rtb221=cneg(cdiv(turf21[ijk1],det));
								rtb222=cdiv(turf11[ijk1],det);
								rdnr11[ip]=cadd(cmul(rtb211,tusf11[ip]),
									cmul(rtb212,tusf21[ip]));
								rdnr12[ip]=cadd(cmul(rtb211,tusf12[ip]),
									cmul(rtb212,tusf22[ip]));
								rdnr21[ip]=cadd(cmul(rtb221,tusf11[ip]),
									cmul(rtb222,tusf21[ip]));
								rdnr22[ip]=cadd(cmul(rtb221,tusf12[ip]),
									cmul(rtb222,tusf22[ip]));
							}
						} else {
			
							/* source elastic,receiver acoustic */
							/* or source and receiver acoustic */
							for (ip=0; ip<block_size; ip++) {
								ijk1=ik1+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);
							}
						}
	
						for (ip=0; ip<block_size; ip++) {
							ijk1=ik1+ip;
							ijk2=ik2+ip;
							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(prs[iz],cadd(cmul(cadd(rtb111,rtb211),
								vuzs1),cmul(cadd(rtb112,rtb212),vuzs2)));
							rtb211=csub(cmul(rurf11[ijk1],p[ip]),cmul(gt[ijk2],
								rurf21[ijk1]));
							rtb212=csub(cmul(rurf12[ijk1],p[ip]),cmul(gt[ijk2],
								rurf22[ijk1]));
							rtb221=cadd(cmul(rurf11[ijk1],gl[ijk2]),
								cmul(rurf21[ijk1],p[ip]));
							rtb222=cadd(cmul(rurf12[ijk1],gl[ijk2]),
								cmul(rurf22[ijk1],p[ip]));
							rtb311=cadd(p[ip],rtb211);
							rtb312=cadd(gt[ijk2],rtb212);
							rtb321=csub(rtb221,gl[ijk2]);
							rtb322=cadd(p[ip],rtb222);
							rtb211=cadd(cmul(rtb311,rtb111),
								cmul(rtb312,rtb121));
							rtb212=cadd(cmul(rtb311,rtb112),
								cmul(rtb312,rtb122));
							rtb221=cadd(cmul(rtb321,rtb111),
								cmul(rtb322,rtb121));
							rtb222=cadd(cmul(rtb321,rtb112),
								cmul(rtb322,rtb122));