psvreflect.c

su 的源代码库

							rvrb111=cdiv(y,det);
							rvrb112=cdiv(rtb112,det);
							rvrb121=cdiv(rtb121,det);
							rvrb122=cdiv(x,det);
							rtb111=cadd(cmul(rvrb111,wrn011),
								cmul(rvrb112,wrn021));
							rtb112=cadd(cmul(rvrb111,wrn012),
								cmul(rvrb112,wrn022));
							rtb121=cadd(cmul(rvrb121,wrn011),
								cmul(rvrb122,wrn021));
							rtb122=cadd(cmul(rvrb121,wrn012),
								cmul(rvrb122,wrn022));
							rtb211=cadd(cmul(ewrd11,rtb111),
								cmul(ewrd12,rtb121));
							rtb212=cadd(cmul(ewrd11,rtb112),
								cmul(ewrd12,rtb122));
							rtb221=cadd(cmul(ewrd21,rtb111),
								cmul(ewrd22,rtb121));
							rtb222=cadd(cmul(ewrd21,rtb112),
								cmul(ewrd22,rtb122));
							rvrb211=cadd(cmplx(1.0,0.0),rtb211);
							rvrb212=rtb212;
							rvrb221=rtb221;
							rvrb222=cadd(cmplx(1.0,0.0),rtb222);
							rtb111=cadd(ewigh11,td0n11[ip]);
							rtb112=cadd(ewigh11,td0n12[ip]);
							rtb121=cadd(ewigh22,td0n21[ip]);
							rtb122=cadd(ewigh22,td0n22[ip]);
							ewd0211=cadd(cmul(rvrb111,rtb111),
								cmul(rvrb112,rtb121));
							ewd0212=cadd(cmul(rvrb111,rtb112),
								cmul(rvrb112,rtb122));
							ewd0221=cadd(cmul(rvrb121,rtb111),
								cmul(rvrb122,rtb121));
							ewd0222=cadd(cmul(rvrb121,rtb112),
								cmul(rvrb122,rtb122));
							td0n11[ip]=cadd(cmul(ewd0211,td11),
								cmul(ewd0221,td12));
							td0n12[ip]=cadd(cmul(ewd0212,td11),
								cmul(ewd0222,td12));
							td0n21[ip]=cadd(cmul(ewd0211,td21),
								cmul(ewd0221,td22));
							td0n22[ip]=cadd(cmul(ewd0212,td21),
								cmul(ewd0222,td22));
							rtb111=cmul(ewigh11,tu11);
							rtb112=cmul(ewigh11,tu12);
							rtb121=cmul(ewigh22,tu21);
							rtb122=cmul(ewigh22,tu22);
							ewtu211=cadd(cmul(rvrb211,rtb111),
								cmul(rvrb212,rtb121));
							ewtu212=cadd(cmul(rvrb211,rtb112),
								cmul(rvrb212,rtb122));
							ewtu221=cadd(cmul(rvrb221,rtb111),
								cmul(rvrb222,rtb121));
							ewtu222=cadd(cmul(rvrb221,rtb112),
								cmul(rvrb222,rtb122));
							tun0p11=cadd(cmul(tun011[ip],ewtu211),
								cmul(tun012[ip],ewtu221));
							tun0p12=cadd(cmul(tun011[ip],ewtu212),
								cmul(tun012[ip],ewtu222));
							tun0p21=cadd(cmul(tun021[ip],ewtu211),
								cmul(tun022[ip],ewtu221));
							tun0p22=cadd(cmul(tun021[ip],ewtu212),
								cmul(tun022[ip],ewtu222));
							rtb111=cadd(cmul(ewrd11,ewd0211),
								cmul(ewd0221,ewrd12));
							rtb112=cadd(cmul(ewrd11,ewd0212),
								cmul(ewd0222,ewrd12));
							rtb121=cadd(cmul(ewrd21,ewd0211),
								cmul(ewd0221,ewrd22));
							rtb122=cadd(cmul(ewrd21,ewd0212),
								cmul(ewd0222,ewrd22));
							rtb211=cadd(cmul(rtb111,tun011[ip]),cmul(rtb121,
								tun012[ip]));
							rtb212=cadd(cmul(rtb112,tun011[ip]),cmul(rtb122,
								tun012[ip]));
							rtb221=cadd(cmul(rtb111,tun021[ip]),cmul(rtb121,
								tun022[ip]));
							rtb222=cadd(cmul(rtb112,tun021[ip]),cmul(rtb122,
								tun022[ip]));
							rd0np11=cadd(rd0n11[ip],rtb211);
							rd0np12=cadd(rd0n12[ip],rtb212);
							rd0np21=cadd(rd0n21[ip],rtb221);
							rd0np22=cadd(rd0n22[ip],rtb222);
							rtb111=cadd(cmul(wrn011,ewtu211),
								cmul(wrn012,ewtu221));
							rtb112=cadd(cmul(wrn011,ewtu212),
								cmul(wrn012,ewtu222));
							rtb121=cadd(cmul(wrn021,ewtu211),
								cmul(wrn022,ewtu221));
							rtb122=cadd(cmul(wrn021,ewtu212),
								cmul(wrn022,ewtu222));
							rtb211=cadd(cmul(td11,rtb111),cmul(td12,rtb121));
							rtb212=cadd(cmul(td11,rtb112),cmul(td12,rtb122));
							rtb221=cadd(cmul(td21,rtb111),cmul(td22,rtb121));
							rtb222=cadd(cmul(td21,rtb112),cmul(td22,rtb122));
							run011[ip]=cadd(ru11,rtb211);
							run012[ip]=cadd(ru12,rtb212);
							run021[ip]=cadd(ru21,rtb221);
							run022[ip]=cadd(ru22,rtb222);
							rd0n11[ip]=rd0np11;
							rd0n12[ip]=rd0np12;
							rd0n21[ip]=rd0np21;
							rd0n22[ip]=rd0np22;
							tun011[ip]=tun0p11;
							tun012[ip]=tun0p12;
							tun021[ip]=tun0p21;
							tun022[ip]=tun0p22;
						}
					} else if (((at1.r!=0.)||(at1.i!=0.))&&((at2.r!=0.)||
						(at2.i!=0.))){
			
						/* case 4: solid-solid */
						if ((al1.r==al2.r)&&(al1.i==al2.i)&&(at1.r==at2.r)
							&&(at1.i==at2.i)) {		
							for(ip=0;ip<block_size;ip++) {
								ijk1=ik1+ip;
								gl1=gl[ijk1];
								gt1=gt[ijk1];
								ewigh11=cexp(cmul(cmplx(0.0,1.0),cmul(t1,cmul(wpie,gl1))));
								ewigh22=cexp(cmul(cmplx(0.0,1.0),cmul(t1,cmul(wpie,gt1))));
								tun0p11=cmul(tun011[ip],ewigh11);
								tun0p12=cmul(tun012[ip],ewigh22);
								tun0p21=cmul(tun021[ip],ewigh11);
								tun0p22=cmul(tun022[ip],ewigh22);
								rtb111=cmul(ewigh11,run011[ip]);
								rtb112=cmul(ewigh11,run012[ip]);
								rtb121=cmul(ewigh22,run021[ip]);
								rtb122=cmul(ewigh22,run022[ip]);
								run011[ip]=cmul(ewigh11,rtb111);
								run012[ip]=cmul(ewigh22,rtb112);
								run021[ip]=cmul(ewigh11,rtb121);
								run022[ip]=cmul(ewigh22,rtb122);
								td0np11=cmul(ewigh11,td0n11[ip]);
								td0np12=cmul(ewigh11,td0n12[ip]);
								td0np21=cmul(ewigh22,td0n21[ip]);
								td0np22=cmul(ewigh22,td0n22[ip]);
								tun011[ip]=tun0p11;
								tun012[ip]=tun0p12;
								tun021[ip]=tun0p21;
								tun022[ip]=tun0p22;
								td0n11[ip]=td0np11;
								td0n12[ip]=td0np12;
								td0n21[ip]=td0np21;
								td0n22[ip]=td0np22;
							}	
						} else {
							d=crmul(csub(cdiv(rho2,cmul(at2,at2)),cdiv(rho1,
								cmul(at1,at1))),2.0);
							for (ip=0;ip<block_size;ip++) {
								ijk1=ik1+ip;
								ijk2=ik2+ip;
								psq=pp[ip];
								gl1=gl[ijk1];
								gt1=gt[ijk1];
								gam1=gam[ijk1];
								gl2=gl[ijk2];
								gt2=gt[ijk2];
								gam2=gam[ijk2];
								ewigh11=cexp(cmul(cmplx(0.0,1.0),cmul(wpie,cmul(gl1,t1))));
								ewigh22=cexp(cmul(cmplx(0.0,1.0),cmul(wpie,cmul(gt1,t1))));
								a=csub(gam2,gam1);
								b=cadd(rho1,a);
								c=csub(rho2,a);
								e=cadd(cmul(gl1,b),cmul(gl2,c));
								f=cadd(cmul(gt1,b),cmul(gt2,c));
								g=csub(a,cmul(gl1,cmul(gt2,d)));
								h=csub(a,cmul(gl2,cmul(gt1,d)));
								z=cdiv(cmplx(1.0,0.0),cadd(cmul(e,f),cmul(g,cmul(psq,h))));
								y=cmul(z,crmul(rho1,2.0));
								td11=cmul(y,cmul(gl1,f));
								td12=cmul(y,cmul(gt1,cmul(p[ip],g)));
								td21=cneg(cmul(y,cmul(gl1,cmul(p[ip],h))));
								td22=cmul(y,cmul(gt1,e));
								y=crmul(cmul(z,cmul(p[ip],cadd(cmul(a,b),
									cmul(gl2,cmul(gt2,cmul(c,d)))))),2.0);
								rd11=cmul(z,csub(cmul(f,csub(cmul(gl1,b),
									cmul(gl2,c))),cmul(cadd(a,cmul(gl1,
									cmul(gt2,d))),cmul(psq,h))));
								rd12=cmul(y,gt1);
								rd21=cneg(cmul(y,gl1));
								rd22=cmul(z,csub(cmul(e,csub(cmul(gt1,b),
									cmul(gt2,c))),cmul(cadd(a,cmul(gt1,
									cmul(gl2,d))),cmul(psq,g))));
								y=cmul(z,crmul(rho2,2.0));
								tu11=cmul(y,cmul(gl2,f));
								tu12=cmul(y,cmul(gt2,cmul(p[ip],h)));
								tu21=cneg(cmul(y,cmul(gl2,cmul(p[ip],g))));
								tu22=cmul(y,cmul(gt2,e));
								y=crmul(cmul(z,cmul(p[ip],cadd(cmul(a,c),
									cmul(gl1,cmul(gt1,cmul(b,d)))))),2.0);
								ru11=cmul(z,csub(cmul(f,csub(cmul(gl2,c),
									cmul(gl1,b))),cmul(cadd(a,cmul(gl2,
									cmul(gt1,d))),cmul(psq,g))));
								ru12=cmul(y,gt2);
								ru21=cneg(cmul(y,gl2));
								ru22=cmul(z,csub(cmul(e,csub(cmul(gt2,c),
									cmul(gt1,b))),cmul(cadd(a,cmul(gl1,
									cmul(gt2,d))),cmul(psq,h))));
								ewrd11=cmul(ewigh11,rd11);
								ewrd12=cmul(ewigh11,rd12);
								ewrd21=cmul(ewigh22,rd21);
								ewrd22=cmul(ewigh22,rd22);
								wrn011=cmul(ewigh11,run011[ip]);
								wrn012=cmul(ewigh11,run012[ip]);
								wrn021=cmul(ewigh22,run021[ip]);
								wrn022=cmul(ewigh22,run022[ip]);
								rtb111=cadd(cmul(wrn011,ewrd11),
									cmul(wrn012,ewrd21));
								rtb112=cadd(cmul(wrn011,ewrd12),
									cmul(wrn012,ewrd22));
								rtb121=cadd(cmul(wrn021,ewrd11),
									cmul(wrn022,ewrd21));
								rtb122=cadd(cmul(wrn021,ewrd12),
									cmul(wrn022,ewrd22));
								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(wrn011,rvrb111),
									cmul(rvrb112,wrn021));
								rtb112=cadd(cmul(wrn012,rvrb111),
									cmul(rvrb112,wrn022));
								rtb121=cadd(cmul(wrn011,rvrb121),
									cmul(rvrb122,wrn021));
								rtb122=cadd(cmul(wrn012,rvrb121),
									cmul(rvrb122,wrn022));
								rtb211=cadd(cmul(ewrd11,rtb111),
									cmul(ewrd12,rtb121));
								rtb212=cadd(cmul(ewrd11,rtb112),
									cmul(ewrd12,rtb122));
								rtb221=cadd(cmul(ewrd21,rtb111),
									cmul(ewrd22,rtb121));
								rtb222=cadd(cmul(ewrd21,rtb112),
									cmul(ewrd22,rtb122));
								rvrb211=cadd(cmplx(1.0,0.0),rtb211);
								rvrb212=rtb212;
								rvrb221=rtb221;
								rvrb222=cadd(cmplx(1.0,0.0),rtb222);
								rtb111=cmul(ewigh11,td0n11[ip]);
								rtb112=cmul(ewigh11,td0n12[ip]);
								rtb121=cmul(ewigh22,td0n21[ip]);
								rtb122=cmul(ewigh22,td0n22[ip]);
								ewd0211=cadd(cmul(rvrb111,rtb111),
									cmul(rvrb112,rtb121));
								ewd0212=cadd(cmul(rvrb111,rtb112),
									cmul(rvrb112,rtb122));
								ewd0221=cadd(cmul(rvrb121,rtb111),
									cmul(rvrb122,rtb121));
								ewd0222=cadd(cmul(rvrb121,rtb112),
									cmul(rvrb122,rtb122));
								td0n11[ip]=cadd(cmul(td11,ewd0211),
									cmul(ewd0221,td12));
								td0n12[ip]=cadd(cmul(td11,ewd0212),
									cmul(ewd0222,td12));
								td0n21[ip]=cadd(cmul(td21,ewd0211),
									cmul(ewd0221,td22));
								td0n22[ip]=cadd(cmul(td21,ewd0212),
									cmul(ewd0222,td22));
								rtb111=cmul(ewigh11,tu11);
								rtb112=cmul(ewigh11,tu12);
								rtb121=cmul(ewigh22,tu21);
								rtb122=cmul(ewigh22,tu22);
								ewtu211=cadd(cmul(rvrb211,rtb111),
									cmul(rvrb212,rtb121));
								ewtu212=cadd(cmul(rvrb211,rtb112),
									cmul(rvrb212,rtb122));
								ewtu221=cadd(cmul(rvrb221,rtb111),
									cmul(rvrb222,rtb121));
								ewtu222=cadd(cmul(rvrb221,rtb112),
									cmul(rvrb222,rtb122));
								tun0p11=cadd(cmul(tun011[ip],ewtu211),
									cmul(tun012[ip],ewtu221));
								tun0p12=cadd(cmul(tun011[ip],ewtu212),	
									cmul(tun012[ip],ewtu222));
								tun0p21=cadd(cmul(tun021[ip],ewtu211),
									cmul(tun022[ip],ewtu221));
								tun0p22=cadd(cmul(tun021[ip],ewtu212),
									cmul(tun022[ip],ewtu222));
								rtb111=cadd(cmul(ewrd11,ewd0211),
									cmul(ewrd12,ewd0221));
								rtb112=cadd(cmul(ewrd11,ewd0212),
									cmul(ewrd12,ewd0222));
								rtb121=cadd(cmul(ewrd21,ewd0211),
									cmul(ewrd22,ewd0221));
								rtb122=cadd(cmul(ewrd21,ewd0212),
									cmul(ewrd22,ewd0222));
								rtb211=cadd(cmul(rtb111,tun011[ip]),cmul(rtb121,
									tun012[ip]));
								rtb212=cadd(cmul(rtb112,tun011[ip]),cmul(rtb122,
									tun012[ip]));
								rtb221=cadd(cmul(rtb111,tun021[ip]),cmul(rtb121,
									tun022[ip]));
								rtb222=cadd(cmul(rtb112,tun021[ip]),cmul(rtb122,
									tun022[ip]));
								rd0np11=cadd(rd0n11[ip],rtb211);
								rd0np12=cadd(rd0n12[ip],rtb212);
								rd0np21=cadd(rd0n21[ip],rtb221);
								rd0np22=cadd(rd0n22[ip],rtb222);
								rtb111=cadd(cmul(wrn011,ewtu211),
									cmul(ewtu221,wrn012));
								rtb112=cadd(cmul(wrn011,ewtu212),
									cmul(ewtu222,wrn012));
								rtb121=cadd(cmul(wrn021,ewtu211),
									cmul(ewtu221,wrn022));
								rtb122=cadd(cmul(wrn021,ewtu212),
									cmul(ewtu222,wrn022));
								rtb211=cadd(cmul(td11,rtb111),
									cmul(td12,rtb121));
								rtb212=cadd(cmul(td11,rtb112),
									cmul(td12,rtb122));
								rtb221=cadd(cmul(td21,rtb111),
									cmul(td22,rtb121));
								rtb222=cadd(cmul(td21,rtb112),
									cmul(td22,rtb122));
								run011[ip]=cadd(ru11,rtb211);
								run012[ip]=cadd(ru12,rtb212);
								run021[ip]=cadd(ru21,rtb221);
								run022[ip]=cadd(ru22,rtb222);
								rd0n11[ip]=rd0np11;
								rd0n12[ip]=rd0np12;
								rd0n21[ip]=rd0np21;
								rd0n22[ip]=rd0np22;
								tun011[ip]=tun0p11;
								tun012[ip]=tun0p12;
								tun021[ip]=tun0p21;
								tun022[ip]=tun0p22;
							}
						}
					}	
				}					/* closes "if" of layers after first one*/

				if (il==lobs[ijk]-2) {
					ik1=ijk*block_size;
					if (acoustic[ijk]==2.||acoustic[ijk]==4.) {
						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];
							tun011[ip]=cmplx(1.0,0.0);
							tun012[ip]=cmplx(0.0,0.0);