psvreflect.c

su 的源代码库

							ruf22=cdiv(ruf22,det);
							ewigh11=cexp(cmul(gl1,cmul(wpie,cmul(cmplx(0.0,1.0),t1))));
							ewigh22=cexp(cmul(gt1,cmul(wpie,cmul(cmplx(0.0,1.0),t1))));
							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);
							td0n11[ip]=ewigh11;
							td0n12[ip]=cmplx(0.0,0.0);
							td0n21[ip]=cmplx(0.0,0.0);
							td0n22[ip]=ewigh22;
							tun011[ip]=ewigh11;
							tun012[ip]=cmplx(0.0,0.0);
							tun021[ip]=cmplx(0.0,0.0);
							tun022[ip]=ewigh22;
							rtb111=cmul(ewigh11,ruf11);	
							rtb112=cmul(ewigh11,ruf12);	
							rtb121=cmul(ewigh22,ruf21);	
							rtb122=cmul(ewigh22,ruf22);	
							run011[ip]=cmul(ewigh11,rtb111);	
							run012[ip]=cmul(ewigh22,rtb112);	
							run021[ip]=cmul(ewigh11,rtb121);	
							run022[ip]=cmul(ewigh22,rtb122);	
						}
					}
				} else {				/* layers after the first one */

					/* case 1: liquid-liquid */
					if ((at1.r==0.)&&(at1.i==0.)&&(at2.r==0.)&&(at2.i==0.)) {
						if ((al1.r==al2.r)&&(al1.i==al2.i)) {
							for (ip=0; ip<block_size; ip++) {
								ijk1=ik1+ip;
								gl1=gl[ijk1];
								ewigh11=cexp(cmul(wpie,cmul(cmplx(0.0,1.0),cmul(t1,gl1))));
								tun0p11=cmul(ewigh11,tun011[ip]);
								tun0p21=cmul(ewigh11,tun021[ip]);
								rtb111=cmul(ewigh11,run011[ip]);
								run011[ip]=cmul(ewigh11,rtb111);
								run012[ip]=cmplx(0.0,0.0);
								run021[ip]=cmplx(0.0,0.0);
								run022[ip]=cmplx(0.0,0.0);
								td0np11=cmul(ewigh11,td0n11[ip]);
								td0np12=cmul(ewigh11,td0n12[ip]);
								tun011[ip]=tun0p11;
								tun012[ip]=cmplx(0.0,0.0);
								tun021[ip]=tun0p21;
								tun022[ip]=cmplx(0.0,0.0);
								td0n11[ip]=td0np11;
								td0n12[ip]=td0np12;
								td0n21[ip]=cmplx(0.0,0.0);
								td0n22[ip]=cmplx(0.0,0.0);
							}
						} else {
							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];
								gam2=gam[ijk2];
								gl2=gl[ijk2];
								gt2=gt[ijk2];
								ewigh11=cexp(cmul(wpie,cmul(cmplx(0.0,1.0),cmul(t1,gl1))));
								a=csub(gam2,gam1);
								b=cadd(rho1,a);
								c=csub(rho2,a);
								e=cadd(cmul(b,gl1),cmul(c,gl2));
								y=cadd(cmul(gl2,rho1),cmul(gl1,rho2));
								td11=cmul(cdiv(gl1,y),crmul(rho1,2.));
								rd11=cdiv(csub(cmul(gl1,rho2),
									cmul(gl2,rho1)),y);
								tu11=cmul(cdiv(gl2,y),crmul(rho2,2.));
								ru11=cneg(rd11);
								ewrd11=cmul(ewigh11,rd11);
								wrn011=cmul(ewigh11,run011[ip]);
								rvrb111=cdiv(cmplx(1.0,0.0),csub(cmplx(1.0,0.0),cmul(wrn011,ewrd11)));
								rvrb211=cadd(cmplx(1.0,0.0),cmul(rvrb111,
									cmul(ewrd11,wrn011)));
								ewd0211=cmul(td0n11[ip],cmul(rvrb111,ewigh11));
								ewd0212=cmul(td0n12[ip],cmul(rvrb111,ewigh11));
								td0n11[ip]=cmul(td11,ewd0211);
								td0n12[ip]=cmul(td11,ewd0212);
								td0n21[ip]=cmplx(0.0,0.0);
								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;
							}
						}
					} else if (((at1.r==0.)&&(at1.i==0.))&&((at2.r!=0.)||
								(at2.i!=0.))) {

						/* case 2: liquid-solid */
						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(wpie,cmul(cmplx(0.0,1.0),cmul(t1,gl1))));
							ewigh22=cmplx(0.0,0.0);
							a=csub(gam2,gam1);
							b=cadd(rho1,a);
							c=csub(rho2,a);
							e=cadd(cmul(gl1,b),cmul(gl2,c));
							d=cdiv(crmul(rho2,2.),cmul(at2,at2));
							f=b;
							g=csub(a,cmul(d,cmul(gl1,gt2)));
							h=cneg(cmul(d,gl2));
							z=cdiv(cmplx(1.0,0.0),cadd(cmul(e,f),cmul(g,cmul(h,psq))));
							y=cmul(cmul(z,gl1),crmul(rho1,2.0));
							td11=cmul(y,f);
							td12=cmplx(0.0,0.0);
							td21=cneg(cmul(y,cmul(h,p[ip])));
							td22=cmplx(0.0,0.0);
							rd11=cmul(z,csub(cmul(f,csub(cmul(gl1,b),
								cmul(gl2,c))),
							cmul(cmul(psq,h),cadd(a,cmul(gl1,cmul(gt2,d))))));
							rd12=cmplx(0.0,0.0);
							rd21=cmplx(0.0,0.0);
							rd22=cmplx(0.0,0.0);
							y=cmul(z,crmul(rho2,2.0));
							tu11=cmul(y,cmul(f,gl2));
							tu12=cmul(p[ip],cmul(y,cmul(gt2,h)));
							tu21=cmplx(0.0,0.0);
							tu22=cmplx(0.0,0.0);
							y=crmul(cmul(p[ip],cmul(gl1,cmul(b,
								cmul(d,z)))),2.0);
							ru11=cmul(z,csub(cmul(f,csub(cmul(gl2,c),
								cmul(gl1,b))),cmul(cmul(gl2,d),cmul(g,psq))));
							ru12=cmul(gt2,y);
							ru21=cneg(cmul(gl2,y));
							ru22=cneg(cmul(z,cadd(cmul(b,e),cmul(cmul(psq,h),
								cadd(a,cmul(gl1,cmul(gt2,d)))))));
							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(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(rtb211,cmplx(1.0,0.0));
							rvrb212=rtb212;
							rvrb221=rtb221;
							rvrb222=cadd(rtb222,cmplx(1.0,0.0));
							rtb111=cmul(ewigh11,td0n11[ip]);
							rtb112=cmul(ewigh11,td0n12[ip]);
							rtb121=cmul(ewigh22,td0n21[ip]);
							rtb122=cmul(ewigh22,td0n22[ip]);
							ewd0211=cadd(cmul(rtb111,rvrb111),
								cmul(rvrb112,rtb121));
							ewd0212=cadd(cmul(rtb112,rvrb111),
								cmul(rvrb112,rtb122));
							ewd0221=cadd(cmul(rtb111,rvrb121),
								cmul(rvrb122,rtb121));
							ewd0222=cadd(cmul(rtb112,rvrb121),
								cmul(rvrb122,rtb122));
							td0n11[ip]=cadd(cmul(td11,ewd0211),
								cmul(td12,ewd0221));
							td0n12[ip]=cadd(cmul(td11,ewd0212),
								cmul(td12,ewd0222));
							td0n21[ip]=cadd(cmul(td21,ewd0211),
								cmul(td22,ewd0221));
							td0n22[ip]=cadd(cmul(td21,ewd0212),
								cmul(td22,ewd0222));
							rtb111=cadd(ewigh11,tu11);
							rtb112=cadd(ewigh11,tu12);
							rtb121=cadd(ewigh22,tu21);
							rtb122=cadd(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(rtb211,rd0n11[ip]);
							rd0np12=cadd(rtb212,rd0n12[ip]);
							rd0np21=cadd(rtb221,rd0n21[ip]);
							rd0np22=cadd(rtb222,rd0n22[ip]);
							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;
						}
					} else if ((at1.r !=0.)&&(at1.i !=0.)&&(at2.r==0.)&&
						(at2.i==0.)) {

						/* case 3: solid-liquid */
						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[ik2];
							gt2=gt[ik2];
							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(b,gl1),cmul(c,gl2));
							d=cdiv(crmul(rho1,-2.0),cmul(at1,at1));
							f=c;
							g=cneg(cmul(gl1,d));
							h=csub(a,cmul(gl2,cmul(d,gt1)));
							z=cdiv(cmplx(1.0,0.0),cadd(cmul(e,f),cmul(g,cmul(h,psq))));
							y=cmul(z,crmul(rho1,2));
							td11=cmul(y,cmul(gl1,f));
							td12=cmul(y,cmul(p[ip],cmul(gt1,g)));
							td21=cmplx(0.0,0.0);
							td22=cmplx(0.0,0.0);
							y=crmul(cmul(z,cmul(c,cmul(d,
								cmul(p[ip],gl2)))),2.0);
							rd11=cmul(z,csub(cmul(f,csub(cmul(gl1,b),
								cmul(gl2,c))),cmul(cmul(psq,h),cmul(gl1,d))));
							rd12=cmul(y,gt1);	
							rd21=cneg(cmul(gl1,y));
							rd22=cneg(cmul(z,cadd(cmul(c,e),cmul(cadd(a,
								cmul(gt1,cmul(gl2,d))),cmul(psq,g)))));
							y=cmul(cmul(z,gl2),crmul(rho2,2.));
							tu11=cmul(y,f);
							tu12=cmplx(0.0,0.0);
							tu21=cneg(cmul(y,cmul(p[ip],y)));
							tu22=cmplx(0.0,0.0);
							ru11=cmul(z,csub(cmul(f,csub(cmul(gl2,c),
								cmul(gl1,b))),cmul(cadd(a,cmul(gl2,
								cmul(d,gt1))),cmul(psq,g))));
							ru12=cmplx(0.0,0.0);
							ru21=cmplx(0.0,0.0);
							ru22=cmplx(0.0,0.0);
							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(ewrd21,wrn012));
							rtb112=cadd(cmul(wrn011,ewrd12),
								cmul(ewrd22,wrn012));
							rtb121=cadd(cmul(wrn021,ewrd11),
								cmul(ewrd21,wrn022));
							rtb122=cadd(cmul(wrn021,ewrd12),
								cmul(ewrd22,wrn022));
							x=csub(cmplx(1.0,0.0),rtb111);
							y=csub(cmplx(1.0,0.0),rtb122);
							det=csub(cmul(x,y),cmul(rtb112,rtb121));