/*EMF.java*/
/**
 * @author Singh T. Junior and Umberto Ravaioli, Maestro Series, EM-Fields
 * @version 1.0
 * Date: May , 1997
 */
//package Maestro.lib.math;
//import Maestro.lib.math.Complex;


public class EMF{
    //Static Constants
    //private static final double epsilon0 = 8.8541878176E-12; //Units: F/m  exact vaue
    private static final double epsilon0 = 8.841941286E-12; //8.8541878176E-12; //Units: F/m  Approximate value
    private static final double mu0 = 1.25663706144E-6; //Units H/m
//Constructor
EMF(){}//Do nothing.

//public static variables
	public static final double beta = 2.0 * Math.PI;

/*-------------------------------------------------------------------------*/
//Given a Gamma, returns normalized Z.
public static final Complex computeZ(Complex Gamma){
	return Complex.Divide(
		Complex.Add(1.0,Gamma),
		Complex.Subtract(1.0,Gamma)
	);
}

/*-------------------------------------------------------------------------*/
//Given a Gamma and Z0, returns Z
public static final Complex computeZ(Complex Gamma, double Z0){
	return Complex.Multiply(computeZ(Gamma),Z0);
}

/*-------------------------------------------------------------------------*/
//Given a Gamma and Complex Z0, returns Z
public static final Complex computeZ(Complex Gamma, Complex Z0){
	return Complex.Multiply(computeZ(Gamma),Z0);
}


/*-------------------------------------------------------------------------*/
//Given normalized Z0, returns Gamma
public static final Complex computeGamma(Complex ZL){
	double Z0 = 1.0;
	return computeGamma(ZL,Z0);
}


/*-------------------------------------------------------------------------*/
//Given ZL and Complex Z0, it returns Gamma
public static final Complex computeGamma(Complex ZL, Complex Z0){
	return Complex.Divide(Complex.Subtract(ZL,Z0),Complex.Add(ZL,Z0));
}


/*-------------------------------------------------------------------------*/
//Given ZL and double Z0, it returns Gamma
public static final Complex computeGamma(Complex ZL, double Z0){
    
	return Complex.Divide(Complex.Subtract(ZL,Z0),Complex.Add(ZL,Z0));
    
}

/*-------------------------------------------------------------------------*/
// for arg = true, X=ZL;
// for arg = false, X=YL
public static final Complex computeGamma(Complex X, boolean arg){
	if(arg){
	    return computeGamma(X,1.0);
	}
	else{
	    return computeGamma(X,1.0,false);
	}
}


/*-------------------------------------------------------------------------*/
// for arg = true, X=ZL;
// for arg = false, X=YL
public static final Complex computeGamma(Complex X, Complex Z0, boolean arg){
	if(arg){
	    
		return Complex.Divide(Complex.Subtract(X,Z0),Complex.Add(X,Z0));
	   
	}
	else{
	    return Complex.Divide(
	       Complex.Subtract(1.0,Complex.Multiply(X,Z0)),
	       Complex.Add(1.0,Complex.Multiply(X,Z0))
	     );
	}
}

/*-------------------------------------------------------------------------*/
// for arg = true, X=ZL;
// for arg = false, X=YL
public static final Complex computeGamma(Complex X, double Z0, boolean arg){
	if(arg){
	
		return Complex.Divide(Complex.Subtract(X,Z0),Complex.Add(X,Z0));
	   
	   
	}
	else{
	     return Complex.Divide(
	       Complex.Subtract(1.0,Complex.Multiply(X,Z0)),
	       Complex.Add(1.0,Complex.Multiply(X,Z0))
	     );
	}
}

/*-------------------------------------------------------------------------*/
//Equation 3.45c, Pozar
//Impedance of a Short stub with length x (units of wavelength) and Z0 = 1
public static final Complex computeZShortStub(double x){
	Complex tmp = new Complex(0.0,Math.tan(2.0*Math.PI*x));
	return tmp;
}

/*-------------------------------------------------------------------------*/
//Equation 3.45c, Pozar
//Impedance of a Short stub with length x (units of wavelength) and double Z0 
public static final Complex computeZShortStub(double x, double Z0){
	Complex tmp =  new Complex(0.0,Z0*Math.tan(2.0*Math.PI*x));
	return tmp;
}

/*-------------------------------------------------------------------------*/
//Equation 3.45c, Pozar, Page 80
//Impedance of a Short stub with length x (units of wavelength) and Complex Z0 
public static final Complex computeZShortStub(double x, Complex Z0){
	Complex tmp =  new Complex(0.0,Math.tan(2.0*Math.PI*x));
	tmp.Multiply(Z0);
	return tmp;
}

/*-------------------------------------------------------------------------*/
//Equation 3.46c, Pozar, Page 80
//Impedance of an Open stub with length x (units of wavelength) and Z0 = 1
public static final Complex computeZOpenStub(double x){
	Complex tmp = new Complex(0.0,-1.0/Math.tan(2.0*Math.PI*x));
	return tmp;
}

/*-------------------------------------------------------------------------*/
//Equation 3.46c, Pozar, Page 80
//Impedance of an Open stub with length x (units of wavelength) and double Z0 
public static final Complex computeZOpenStub(double x, double Z0){
	Complex tmp = new Complex(0.0,-Z0/Math.tan(2.0*Math.PI*x));
	return tmp;
}

/*-------------------------------------------------------------------------*/
//Equation 3.46c, Pozar, Page 80
//Impedance of an Open stub with length x (units of wavelength) and Complex Z0 
public static final Complex computeZOpenStub(double x, Complex Z0){
	Complex tmp = new Complex(0.0,-1.0/Math.tan(2.0*Math.PI*x));
	tmp.Multiply(Z0);
	return tmp;
}


/*-------------------------------------------------------------------------*/
/**
 * Compute Gamma at certain distance from Load
 * @param Complex GammaL (reflection coef at load), double beta (wave vector), double x (distance)
 * @return Complex GammaL(x)
 * David Pozar, Microwave Engineering, Eq.3.42, P.79
 */
//Loss Less
public static final Complex computeGammaAt(Complex GammaL, double x){
	double farg = 2*beta*x;
	return Complex.Multiply(new Complex(Math.cos(farg),-Math.sin(farg)),GammaL);
}
/*-------------------------------------------------------------------------*/
//Lossy
public static final Complex computeGammaAt(Complex GammaL, double alpha, double x){
	return Complex.Multiply(computeGammaAt(GammaL,x),Math.exp(-2*alpha*x));
}
/*-------------------------------------------------------------------------*/
/**
 * Compute Zin at a distance x from the load. GammaL is load's Gamma.
 * David Pozar, Microwave Engineering, Eq.3.43, P.79
 */
private static final Complex computeZinAt(Complex GammaL, double Z0, double x){
	Complex GammaShort = new Complex(-1.0,0.0);
	Complex GammaOpen = new Complex(1.0,0.0);
	Complex Small = new Complex(0.0,0.0);
	Complex Large = new Complex(0.0,-1.0E135);
	
	double farg = 2*beta*x;
	
	double farg2 = beta*x;
	double fcos = Math.cos(farg2);
	double fsin = Math.sin(farg2);
	double ftan = fsin/fcos;
	double fcotan = fcos/fsin;
	// FIX for short and open circuits - 7/31/97 - Umberto
	if(GammaL.equals(GammaShort)){
	Complex tmp = new Complex(0.0,ftan);
	Complex refl = computeGammaAt(GammaL,x);
	double phase = Complex.Arg2(refl,1); 
	double dif_phase = Math.abs(phase) - 180.0;
	
	
	    if(Math.abs(phase)< 0.00001){return Large;}
	    
	    if(Math.abs(dif_phase)<0.00001){return Small;}
	    
	    else{
	    return Complex.Multiply(tmp,Z0);}
	}
	
	if(GammaL.equals(GammaOpen)){
	Complex tmp = new Complex(0.0,-fcotan);
	Complex refl = computeGammaAt(GammaL,x);
	double phase = Complex.Arg2(refl,1); 
	double dif_phase = Math.abs(phase) - 180.0;
	
	
	    if(Math.abs(phase)< 0.00001){
		return Large;}
	    
	    if(Math.abs(dif_phase)<0.00001){
		return Small;}
	    
	    else{
		return Complex.Multiply(tmp,Z0);}
	}
	
	// FIX ends here
	
	else{
	Complex tmp = Complex.Multiply(new Complex(Math.cos(farg),-Math.sin(farg)),GammaL);
	return Complex.Multiply(Complex.Divide(Complex.Add(1.0,tmp),Complex.Subtract(1.0,tmp)),Z0);
	}
}

/*-------------------------------------------------------------------------*/
// Computes Zin at a distance x from load. Z0 is complex in this case. GammaL is load's Gamma.
private static final Complex computeZinAt(Complex GammaL, Complex Z0, double x){
	double farg = 2*beta*x;
	Complex tmp = Complex.Multiply(new Complex(Math.cos(farg),-Math.sin(farg)),GammaL);
	return Complex.Multiply(Complex.Divide(Complex.Add(1.0,tmp),Complex.Subtract(1.0,tmp)),Z0);
}

/*-------------------------------------------------------------------------*/
/**
 * Compute Zin
 * If method = true, use GammaL to compute Zin
 * Example: computeZinAt(GammaL,Z0,x,true);
 * If method = false, use ZL to compute Zin
 * Example: computeZinAt(ZL,Z0,x,false);
 */
public static final Complex computeZinAt(Complex jtemp, double Z0, double x, boolean method){
	if(method){
		//jtemp is GammaL, method=true
		return computeZinAt(jtemp,Z0,x);
	}
	else{
		//jtemp is ZL, method=false;
		return computeZinAt(computeGamma(jtemp,Z0),Z0,x);
	}
}

/*-------------------------------------------------------------------------*/
// If jtemp is GammaL, then method=true;
// If jtemp is ZL, then method=false;
public static final Complex computeZinAt(Complex jtemp, Complex Z0, double x, boolean method){
	if(method){
		//jtemp is GammaL, method=true
		return computeZinAt(jtemp,Z0,x);
	}
	else{
		//jtemp is ZL, method=false;
		return computeZinAt(computeGamma(jtemp,Z0,true),Z0,x);
	}
}


/*-------------------------------------------------------------------------*/
/**
    Returns Yin at a distance x from the load with 3 stubs in the line.

    If jtemp is GammaL, then method = true;
    If jtemp is YL, then method = false;
*/


/*-------------------------------------------------------------------------*/
/**
 * Compute Yin at certain distance from GammaL
 * David Pozar, Microwave Engineering, Eq.3.43, P.79
 */
private static final Complex computeYinAt(Complex GammaL, double Z0, double x){
	Complex GammaShort = new Complex(-1.0,0.0);
	Complex GammaOpen = new Complex(1.0,0.0);
	Complex Small = new Complex(0.0,0.0);
	Complex Large = new Complex(0.0,-1E135);
	
	double farg = 2*beta*x;
	
	double farg2 = beta*x;
	double fcos = Math.cos(farg2);
	double fsin = Math.sin(farg2);
	double ftan = fsin/fcos;
	double fcotan = fcos/fsin;
	// FIX for short and open circuits - 7/31/97 - Umberto
	if(GammaL.equals(GammaShort)){
	Complex tmp = new Complex(0.0,-fcotan);
	Complex refl = computeGammaAt(GammaL,x);
	double phase = Complex.Arg2(refl,1); 
	double dif_phase = Math.abs(phase) - 180.0;
	
	
	    if(Math.abs(phase)< 0.00001){
		return Small;}
	    
	    if(Math.abs(dif_phase)<0.00001){
		return Large;}
	    
	    else{
		return Complex.Divide(tmp,Z0);}
	}
	
	if(GammaL.equals(GammaOpen)){
	Complex tmp = new Complex(0.0,ftan);
	Complex refl = computeGammaAt(GammaL,x);
	double phase = Complex.Arg2(refl,1); 
	double dif_phase = Math.abs(phase) - 180.0;
	
	
	    if(Math.abs(phase)< 0.00001){
		return Small;}
	    
	    if(Math.abs(dif_phase)<0.00001){
		return Large;}
	    
	    else{
		return Complex.Divide(tmp,Z0);}
	}
	
	
	// FIX ends here
	
    	else{
	    Complex tmp = Complex.Multiply(new Complex(Math.cos(farg),-Math.sin(farg)),GammaL);
	    return Complex.Divide(Complex.Divide(Complex.Subtract(1.0,tmp),Complex.Add(1.0,tmp)),Z0);
	}
}

/*-------------------------------------------------------------------------*/
// Compute Yiin at certain distance x from GammaL, with Complex Z0.
private static final Complex computeYinAt(Complex GammaL, Complex Z0, double x){
	double farg = 2*beta*x;
	Complex tmp = Complex.Multiply(new Complex(Math.cos(farg),-Math.sin(farg)),GammaL);
	return Complex.Divide(Complex.Divide(Complex.Subtract(1.0,tmp),Complex.Add(1.0,tmp)),Z0);
}

/*-------------------------------------------------------------------------*/
/**
 * Compute Yin
 * If method = true, use GammaL to compute Yin
 * Example: computeYinAt(GammaL,Z0,x,true);
 * If method = false, use YL to compute Zin
 * Example: computeYinAt(YL,Z0,x,false);
 */
public static final Complex computeYinAt(Complex jtemp, double Z0, double x, boolean method){
	if(method){
		//jtemp is GammaL, method=true
		return computeYinAt(jtemp,Z0,x);
	}
	else{
		//jtemp is YL, method=false;
		return computeYinAt(computeGamma(jtemp,Z0,false),Z0,x);
	}
}


/*-------------------------------------------------------------------------*/
public static final Complex computeYinAt(Complex jtemp, Complex Z0, double x, boolean method){
	if(method){
		//jtemp is GammaL, method=true
		return computeYinAt(jtemp,Z0,x);
	}
	else{
		//jtemp is YL, method=false;
		return computeYinAt(computeGamma(jtemp,Z0),Z0,x);
	}
}

/*-------------------------------------------------------------------------*/
/**
 * Equation 8.73, Rao
 */
public static final Complex computeVat(Complex ZL, double Z0, Complex VPlus, double x){
	Complex GammaL=computeGamma(ZL,Z0);
	double  farg = beta*x;
	Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	Complex tmp2 = Complex.Conjugate(tmp1);
	tmp2.Multiply(GammaL);
	return Complex.Multiply(Complex.Add(tmp1,tmp2),VPlus);
}

/*-------------------------------------------------------------------------*/
/**
 * If method = true, X = GammaL
 * If method = false, X = ZL
 */
public static final Complex computeVat(Complex Xtemp, double Z0, Complex VPlus, double x, boolean method){
	if(method){
	    double  farg = beta*x;
	    Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	    Complex tmp2 = Complex.Multiply(Complex.Conjugate(tmp1),Xtemp);
	    return Complex.Multiply(Complex.Add(tmp1,tmp2),VPlus);
	}
	else{
	    return computeVat(Xtemp,Z0,VPlus,x);
	}
}

public static final Complex computeVincidentat(Complex ZL, double Z0, Complex VPlus, double x){
	Complex GammaL=computeGamma(ZL,Z0);
	double  farg = beta*x;
	Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	
	return Complex.Multiply(tmp1,VPlus);
}

public static final Complex computeVreflectedat(Complex ZL, double Z0, Complex VPlus, double x){
	Complex GammaL=computeGamma(ZL,Z0);
	double  farg = beta*x;
	Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	Complex tmp2 = Complex.Conjugate(tmp1);
	tmp2.Multiply(GammaL);
	return Complex.Multiply(tmp2,VPlus);
}
/*-------------------------------------------------------------------------*/
/**
 * If method = true, X = GammaL
 * If method = false, X = ZL
 */

public static final Complex computeIat(Complex Xtemp, double Z0, Complex VPlus, double x, boolean method){
	if(method){
	    double  farg = beta*x;
	    Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	    Complex tmp2 = Complex.Multiply(Complex.Conjugate(tmp1),Xtemp);
	    return Complex.Divide(Complex.Multiply(Complex.Subtract(tmp1,tmp2),VPlus),Z0);
	}
	else{
	    return computeIat(Xtemp,Z0,VPlus,x);
	}
}

/*-------------------------------------------------------------------------*/
public static final Complex computeIat(Complex ZL, double Z0, Complex VPlus, double x){
	Complex GammaL=computeGamma(ZL,Z0);
	double  farg = beta*x;
	Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	Complex tmp2 = Complex.Conjugate(tmp1);
	tmp2.Multiply(GammaL);
	return Complex.Divide(Complex.Multiply(Complex.Subtract(tmp1,tmp2),VPlus),Z0);
}

public static final Complex computeIincidentat(Complex ZL, double Z0, Complex VPlus, double x){
	Complex GammaL=computeGamma(ZL,Z0);
	double  farg = beta*x;
	Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	
	return Complex.Divide(Complex.Multiply(tmp1,VPlus),Z0);
}

public static final Complex computeIreflectedat(Complex ZL, double Z0, Complex VPlus, double x){
	Complex GammaL=computeGamma(ZL,Z0);
	double  farg = beta*x;
	Complex tmp1 = new Complex(0.0,0.0);
        Complex tmp2 = new Complex(Math.cos(farg),-Math.sin(farg));
	//Complex tmp2 = Complex.Conjugate(tmp1);
	tmp2.Multiply(GammaL);
	return Complex.Divide(Complex.Multiply(Complex.Subtract(tmp1,tmp2),VPlus),Z0);
}
/*-------------------------------------------------------------------------*/
/**
 * Equation 8.73, Rao
 */
public static final Complex computeVat(Complex ZL, Complex Z0, Complex VPlus, double x){
	Complex GammaL=computeGamma(ZL,Z0);
	double  farg = beta*x;
	Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	Complex tmp2 = Complex.Conjugate(tmp1);
	tmp2.Multiply(GammaL);
	return Complex.Multiply(Complex.Add(tmp1,tmp2),VPlus);
}

/*-------------------------------------------------------------------------*/
/**
 * If method = true, X = GammaL
 * If method = false, X = ZL
 */
public static final Complex computeVat(Complex Xtemp, Complex Z0, Complex VPlus, double x, boolean method){
	if(method){
	    double  farg = beta*x;
	    Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	    Complex tmp2 = Complex.Multiply(Complex.Conjugate(tmp1),Xtemp);
	    return Complex.Multiply(Complex.Add(tmp1,tmp2),VPlus);
	}
	else{
	    return computeVat(Xtemp,Z0,VPlus,x);
	}
}

/*-------------------------------------------------------------------------*/
/**
 * If method = true, X = GammaL
 * If method = false, X = ZL
 */

public static final Complex computeIat(Complex Xtemp, Complex Z0, Complex VPlus, double x, boolean method){
	if(method){
	    double  farg = beta*x;
	    Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	    Complex tmp2 = Complex.Multiply(Complex.Conjugate(tmp1),Xtemp);
	    return Complex.Divide(Complex.Multiply(Complex.Subtract(tmp1,tmp2),VPlus),Z0);
	}
	else{
	    return computeIat(Xtemp,Z0,VPlus,x);
	}
}

/*-------------------------------------------------------------------------*/
public static final Complex computeIat(Complex ZL, Complex Z0, Complex VPlus, double x){
	Complex GammaL=computeGamma(ZL,Z0);
	double  farg = beta*x;
	Complex tmp1 = new Complex(Math.cos(farg),Math.sin(farg));
	Complex tmp2 = Complex.Conjugate(tmp1);
	tmp2.Multiply(GammaL);
	return Complex.Divide(Complex.Multiply(Complex.Subtract(tmp1,tmp2),VPlus),Z0);
}

/*-------------------------------------------------------------------------*/
//Returns Voltage at the generator terminals for a case without stubs
public static final Complex computeVPlus(Generator g, Complex ZL, double Z0, double lineLength){
	Complex Zin;
	Complex one = new Complex(1.0,0.0);
	Zin = computeZinAt(ZL,Z0,lineLength,false);
	if(Zin.Magnitude2()>0.0){
	  return Complex.Divide(
	    Complex.Divide(
		    Complex.Multiply(g.getVg(),Zin),
		    Complex.Add(g.getZg(),Zin)
	    ),EMF.computeVat(ZL,Z0,one,lineLength));
	}
	else{
	   return Complex.Divide(
	    Complex.Divide(g.getVg(),g.getZg()),EMF.computeIat(ZL,Z0,one,lineLength));	
	}
}

/*-------------------------------------------------------------------------*/
//Returns Voltage at the generator terminals for a case without stubs
public static final Complex computeVPlus(Generator g, Complex ZL, Complex Z0, double lineLength){
	Complex Zin;
	Complex one = new Complex(1.0,0.0);
	Zin = computeZinAt(ZL,Z0,lineLength,false);
	if(Zin.Magnitude2()>0.0){
	  return Complex.Divide(
	    Complex.Divide(
		    Complex.Multiply(g.getVg(),Zin),
		    Complex.Add(g.getZg(),Zin)
	    ),EMF.computeVat(ZL,Z0,one,lineLength));
	}
	else{
	   return Complex.Divide(
	    Complex.Divide(g.getVg(),g.getZg()),EMF.computeIat(ZL,Z0,one,lineLength));	
	}
}



/*-------------------------------------------------------------------------*/
//Returns Voltage at the input terminals, given generator and input impedance at the terminal
public static final Complex voltageDivider(Generator g, Complex Zx){
    	return Complex.Divide(
		    Complex.Multiply(g.getVg(),Zx),
			Complex.Add(Zx,g.getZg())
	);
}

/*-------------------------------------------------------------------------*/
//Returns Voltage at the terminals, given Vg, Zg (for generator), and Zx is input impedance at the terminal
public static final Complex voltageDivider(Complex Vg, Complex Zg, Complex Zx){
    	return Complex.Divide(
		    Complex.Multiply(Vg,Zx),
			Complex.Add(Zx,Zg)
	);
}

/*-------------------------------------------------------------------------*/
public static final Complex Inv(Complex Z){
	Complex Y;
	double mag2;
	
	if(Z.Imaginary()==0.0 && Z.Real() != 0.0){
	    Y = new Complex(1.0/Z.Real(),0.0);
	    return Y;
	}
	if(Z.Real()==0.0 && Z.Imaginary() != 0.0){
	    Y = new Complex(0.0,-1.0/Z.Imaginary());
	    return Y;
	}
	mag2 = Z.Magnitude2();	
	if(mag2>0.0){
	    Y = new Complex(Z.Real()/mag2,-Z.Imaginary()/mag2);
	}
	//else {Y=new Complex(0.0,Double.NEGATIVE_INFINITY); }
	else { Y = new Complex(0.0,-1.0E135); }
	return Y;
}
/*-------------------------------------------------------------------------*/
public static final Complex Series(Complex Z1, Complex Z2){
	return Complex.Add(Z1,Z2);
}	
/*-------------------------------------------------------------------------*/
public static final Complex Parallel(Complex Z1, Complex Z2){
	return Inv(Complex.Add(Inv(Z1),Inv(Z2)));
}

	

/*-------------------------------------------------------------------------*/
public static final double computeSWR(Complex Gamma){
	return (1.0+Gamma.Magnitude())/(1.0-Gamma.Magnitude());
}
/*-------------------------------------------------------------------------*/
/**
 * Microwave Engineering, David M. Pozar, Equation 3.72, Page 98.
 *
 */
public static final double PowerDelivered(Generator G, Complex ZL, double Z0, double length){
	double PDel;
	Complex Zin;
	Zin = computeZinAt(ZL,Z0,length,false);
	PDel =  0.5*Complex.Magnitude2(G.getVg())*Zin.Real()/
	    (Math.pow(Zin.Real()+Complex.Real(G.getZg()),2.0)+
	     Math.pow(Zin.Imaginary()+Complex.Imaginary(G.getZg()),2.0)
            );
	return PDel;
}
/*-------------------------------------------------------------------------*/
/**
 * Microwave Engineering, David M. Pozar, Equation 3.90, 3.93a, 3.93b, page 107.
 *
 */
public static final Complex computePropagator(double R, double L, double G, double C, double w, boolean isLowLoss){
	Complex gamma;
	if(!isLowLoss){
		gamma=Complex.Pow(Complex.Multiply(new Complex(R,w*L),new Complex(G,w*C)),0.5);
		
	}
	else{
		gamma=new Complex(0.5*(R*Math.sqrt(C/L)+G*Math.sqrt(L/C)),
				w*Math.sqrt(L*C));
	}
	return gamma;
}
/*-------------------------------------------------------------------------*/
/**
 * Microwave Engineering, David M. Pozar, Equation 3.100, 3.101, and 3.102, page 110.
 */
public static final double[] PowerDelivered(Generator G,Complex VPlus,
                                  Complex ZL,double alpha, double Z0, double length){
	double P[] = new double[7];
	Complex GammaL;
	Complex Zin;
	Complex Ztot;
	Complex Iin;
	GammaL=computeGamma(ZL,Z0);
	Zin = computeZinAt(ZL,Z0,length,false);
	Ztot = Complex.Add(Zin,G.getZg());
	Iin = Complex.Divide(G.getVg(),Ztot);
	//Time-average (Absorbed) Power at the input of the transmission line.
	P[0] = (VPlus.Magnitude2()/(2.0*Z0))*
	      (Math.exp(2.0*alpha*length)-GammaL.Magnitude2()*Math.exp(-2.0*alpha*length));
	//Time-average (Absorbed) Power delivered to the load
	P[1] = (VPlus.Magnitude2()/(2.0*Z0))*
	       (1.0-GammaL.Magnitude2());
	//Power loss in the line
	P[2] = P[0] - P[1];
	//Time-Average Power Injected into the Transmission Line by the incoming Wave
	P[3] = (VPlus.Magnitude2()/(2.0*Z0))*Math.exp(2.0*alpha*length);
	//Time-Average Power Reflected by the load
	P[4] = (VPlus.Magnitude2()/(2.0*Z0))*GammaL.Magnitude2()*Math.exp(-2.0*alpha*length);
	//Time-Average Power Absorbed by Generator Impedance
	P[5] = 0.5*Complex.Real(Complex.Multiply(Complex.Multiply(G.getZg(),Iin),Complex.Conjugate(Iin)));
	//Total Time-Average Power Spent by the Generator
	P[6] = 0.5*Complex.Real(Complex.Multiply(G.getVg(),Complex.Conjugate(Iin)));
    	return P;
} 
/*-------------------------------------------------------------------------*/
public static final double TimeAveragedPower(Complex V, Complex I){
    return 0.5*Complex.Real(Complex.Multiply(V,Complex.Conjugate(I)));
}
/*-------------------------------------------------------------------------*/
public static final Complex Power(Complex V, Complex I){
    return Complex.Multiply(V,Complex.Conjugate(I));
}
/*-------------------------------------------------------------------------*/
public static final double getPhaseVelocity(double epsilon_r, double mu_r, double conductivity,
    double angular_frequency){
/* Reference:
 * Elements of Engineering Electromagnetics, Fourth Edition
 * Nannapaneni Narayana Rao, Prentice Hall, ISBN 0-13-948746-8
 */
 //Equation 6.67, Page 326
    return Math.sqrt(2.0/(mu0*mu_r*epsilon_r*epsilon0))*
	 Math.pow(
	    Math.pow(
	       1.0+Math.pow((conductivity/(angular_frequency*epsilon_r*epsilon0)),2.0)
	       ,0.5)+1.0
	,-0.5);
}
/*-------------------------------------------------------------------------*/
public static final double getWaveLength(double epsilon_r, double mu_r, double conductivity,
    double angular_frequency){
    //N.N. Rao, Equation 6.68, Page 326
    return (2.0*Math.PI/angular_frequency)*
	     getPhaseVelocity(epsilon_r,mu_r,conductivity,angular_frequency);
}

/*-------------------------------------------------------------------------*/
public static final double getAlpha(double epsilon_r, double mu_r, double conductivity, 
    double angular_frequency){
    //N.N. Rao, Equation 6.65, Page 324
     return  angular_frequency * 
	 Math.sqrt(mu0*mu_r*epsilon_r*epsilon0/2.0)*
	 Math.pow(
	    Math.pow(
	       1.0+Math.pow((conductivity/(angular_frequency*epsilon_r*epsilon0)),2.0)
	       ,0.5)-1.0
	,0.5);
}
/*-------------------------------------------------------------------------*/
public static final double getBeta(double epsilon_r, double mu_r, double conductivity, 
    double angular_frequency){
    //N.N. Rao, Equation 6.65, Page 324
     return angular_frequency * 
	 Math.sqrt(mu0*mu_r*epsilon_r*epsilon0/2.0)*
	 Math.pow(
	    Math.pow(
	       1.0+Math.pow((conductivity/(angular_frequency*epsilon_r*epsilon0)),2.0)
	       ,0.5)+1.0
	,0.5);
}
/*-------------------------------------------------------------------------*/
public static final Complex getWaveImpedance(double epsilon_r, double mu_r, double conductivity, 
    double angular_frequency){
    Complex tmp1, tmp2;
    tmp1 = new Complex(0.0,angular_frequency*mu_r*mu0);
    tmp2 = new Complex(conductivity,angular_frequency*epsilon_r*epsilon0);
    return Complex.Pow(Complex.Divide(tmp1,tmp2),0.5);
}
/*-------------------------------------------------------------------------*/
public static final  Complex getOblique_Reflection_Coefficient(double theta1, double epsilon_r1, double epsilon_r2,
     double mu_r1, double mu_r2, double conductivity1, double conductivity2, double angular_frequency,
     boolean parallel){
	Complex tmp, costheta2;
	double beta1, beta2, theta2,factor1, factor2;
	Complex wave_impedance1, wave_impedance2;
	wave_impedance1 = getWaveImpedance(epsilon_r1,mu_r1,conductivity1,angular_frequency);
	wave_impedance2 = getWaveImpedance(epsilon_r2,mu_r2,conductivity2,angular_frequency);
	beta1 = getBeta(epsilon_r1,mu_r1,conductivity1,angular_frequency);
	beta2 = getBeta(epsilon_r2,mu_r2,conductivity2,angular_frequency);
	//factor1 = Complex.Magnitude(wave_impedance2)*mu_r1/Complex.Magnitude(wave_impedance1)/mu_r2;
	factor1 = Math.sqrt(mu_r1*epsilon_r1/(mu_r2*epsilon_r2));
	
	//theta2 = Math.asin(beta1*Math.sin(theta1)/beta2);
	theta2 = Math.asin(Math.sin(theta1)*factor1);
	if(parallel){
	    tmp = Complex.Divide(
		Complex.Subtract(
		    Complex.Multiply(wave_impedance2,Math.cos(theta2)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta1))
		),
		Complex.Add(
		    Complex.Multiply(wave_impedance2,Math.cos(theta2)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta1))
		)
	    );
	}
	else {
	    tmp = Complex.Divide(
		Complex.Subtract(
		    Complex.Multiply(wave_impedance2,Math.cos(theta1)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta2))
		),
		Complex.Add(
		    Complex.Multiply(wave_impedance2,Math.cos(theta1)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta2))
		)
	    );
	}
	if((mu_r1*epsilon_r1) <= (mu_r2*epsilon_r2)){
		    return tmp;
	  }
	  else{
		    factor2 = Math.sin(theta1)*factor1;
		    if(factor2 > 1.0){
		    
			costheta2 = new Complex(0.0,Math.sqrt(Math.sin(theta1)*Math.sin(theta1)*factor1*factor1-1.0));
			
			if(parallel){
			       tmp = Complex.Divide(
				    Complex.Subtract(
				    Complex.Multiply(wave_impedance2,costheta2),
				    Complex.Multiply(wave_impedance1,Math.cos(theta1))
				    ),
				    Complex.Add(
				    Complex.Multiply(wave_impedance2,costheta2),
				    Complex.Multiply(wave_impedance1,Math.cos(theta1))
				    )
				);
			}
			else {
				    tmp = Complex.Divide(
					Complex.Subtract(
					Complex.Multiply(wave_impedance2,Math.cos(theta1)),
					Complex.Multiply(wave_impedance1,costheta2)
					),
					Complex.Add(
					Complex.Multiply(wave_impedance2,Math.cos(theta1)),
					Complex.Multiply(wave_impedance1,costheta2)
					)
				    );
			}
			//tmp = new Complex(1.0,0.0);
			return tmp;
		    }
		    else{   
			return tmp;
		    }
	  }
    }
    /*-------------------------------------------------------------------------*/
    public static final  Complex getOblique_Transmission_Coefficient(double theta1, double epsilon_r1, double epsilon_r2,
     double mu_r1, double mu_r2, double conductivity1, double conductivity2, double angular_frequency,
     boolean parallel){
	Complex tmp, costheta2;
	double beta1, beta2, theta2, factor1, factor2;
	Complex wave_impedance1, wave_impedance2;
	wave_impedance1 = getWaveImpedance(epsilon_r1,mu_r1,conductivity1,angular_frequency);
	wave_impedance2 = getWaveImpedance(epsilon_r2,mu_r2,conductivity2,angular_frequency);
	beta1 = getBeta(epsilon_r1,mu_r1,conductivity1,angular_frequency);
	beta2 = getBeta(epsilon_r2,mu_r2,conductivity2,angular_frequency);
	//factor1 = Complex.Magnitude(wave_impedance2)*mu_r1/Complex.Magnitude(wave_impedance1)/mu_r2;
	factor1 = Math.sqrt(mu_r1*epsilon_r1/(mu_r2*epsilon_r2));
	
	//theta2 = Math.asin(beta1*Math.sin(theta1)/beta2);
	theta2 = Math.asin(Math.sin(theta1)*factor1);
    	if(parallel){
	    tmp = Complex.Divide(
		Complex.Multiply(wave_impedance2,2*Math.cos(theta1)),
		Complex.Add(
		    Complex.Multiply(wave_impedance2,Math.cos(theta2)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta1))
		)
	    );
	}
	else {
	    tmp = Complex.Divide(
		Complex.Multiply(wave_impedance2,2.0*Math.cos(theta1)),
		Complex.Add(
		    Complex.Multiply(wave_impedance2,Math.cos(theta1)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta2))
		)
	    );
	}
	if((mu_r1*epsilon_r1) <= (mu_r2*epsilon_r2)){
		    return tmp;
	  }
	  else{
		    //factor2 = Math.sin(theta1)*Math.sin(theta1)*factor1*factor1;
		    factor2 = Math.sin(theta1)*factor1;
		    if(factor2 > 1.0){
		    
			costheta2 = new Complex(0.0,Math.sqrt(Math.sin(theta1)*Math.sin(theta1)*factor1*factor1-1.0));
			
			if(parallel){
				    tmp = Complex.Divide(
					Complex.Multiply(wave_impedance2,2*Math.cos(theta1)),
					Complex.Add(
					Complex.Multiply(wave_impedance2,costheta2),
					Complex.Multiply(wave_impedance1,Math.cos(theta1))
					)
				    );
			}
			else {
				    tmp = Complex.Divide(
					Complex.Multiply(wave_impedance2,2.0*Math.cos(theta1)),
					Complex.Add(
					Complex.Multiply(wave_impedance2,Math.cos(theta1)),
					Complex.Multiply(wave_impedance1,costheta2)
					)
				    );
			}
		    
		    
			//tmp = new Complex(0.0,0.0);
			return tmp;
		    }
		    else{   
			return tmp;
		    }
	  }
    }
    /*-------------------------------------------------------------------------*/
    public static final void getOblique_Reflection_Coefficient(double epsilon_r1, double epsilon_r2,
     double mu_r1, double mu_r2, double conductivity1, double conductivity2, double angular_frequency,
     boolean parallel, double array[], double x[]){
	Complex tmp;
	double beta1, beta2, theta1, theta2, factor1, factor2;
	Complex wave_impedance1, wave_impedance2;
	wave_impedance1 = getWaveImpedance(epsilon_r1,mu_r1,conductivity1,angular_frequency);
	wave_impedance2 = getWaveImpedance(epsilon_r2,mu_r2,conductivity2,angular_frequency);
	beta1 = getBeta(epsilon_r1,mu_r1,conductivity1,angular_frequency);
	beta2 = getBeta(epsilon_r2,mu_r2,conductivity2,angular_frequency);
	//factor1 = Complex.Magnitude(wave_impedance2)*mu_r1/Complex.Magnitude(wave_impedance1)/mu_r2;
	factor1 = Math.sqrt(mu_r1*epsilon_r1/(mu_r2*epsilon_r2));
	
	for(int i =  0; i <x.length; i++){
	  theta1 = i * Math.PI/180.0;
	  //MARKER THETA2
	  //theta2 = Math.asin(beta1*Math.sin(theta1)/beta2);
	  theta2 = Math.asin(Math.sin(theta1)*factor1);
	  
	  //System.out.println("theta1 = "+theta1+"  theta2 = "+theta2+"   eta1 = "+wave_impedance1+"   eta2 = "+wave_impedance2);
	  x[i] = theta1;
	  if(parallel){
	    tmp = Complex.Divide(
		Complex.Subtract(
		    Complex.Multiply(wave_impedance2,Math.cos(theta2)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta1))
		),
		Complex.Add(
		    Complex.Multiply(wave_impedance2,Math.cos(theta2)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta1))
		)
	    );
	  }
	  else {
	    tmp = Complex.Divide(
		Complex.Subtract(
		    Complex.Multiply(wave_impedance2,Math.cos(theta1)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta2))
		),
		Complex.Add(
		    Complex.Multiply(wave_impedance2,Math.cos(theta1)),
		    Complex.Multiply(wave_impedance1,Math.cos(theta2))
		)
	    );
	  }
	  if((mu_r1*epsilon_r1) <= (mu_r2*epsilon_r2)){
		    array[i] = tmp.Magnitude()*tmp.Magnitude();
	  }
	  else{
		    //factor2 = Math.sin(theta1)*Math.sin(theta1)*factor1*factor1;
		    factor2 = Math.sin(theta1)*factor1;
		    if(factor2 > 1.0){
			array[i] = 1.0;
		    }
		    else{   
			array[i] = tmp.Magnitude()*tmp.Magnitude();
		    }
	  }
      }//end of loop in i;
    }
    /*-------------------------------------------------------------------------*/
    public static final  void getOblique_Transmission_Coefficient(double epsilon_r1, double epsilon_r2,
     double mu_r1, double mu_r2, double conductivity1, double conductivity2, double angular_frequency,
     boolean parallel, double array[], double x[]){
	Complex tmp;
	double beta1, beta2, theta1, theta2, factor1, factor2, tmp2;
	Complex wave_impedance1, wave_impedance2;
	wave_impedance1 = getWaveImpedance(epsilon_r1,mu_r1,conductivity1,angular_frequency);
	wave_impedance2 = getWaveImpedance(epsilon_r2,mu_r2,conductivity2,angular_frequency);
	beta1 = getBeta(epsilon_r1,mu_r1,conductivity1,angular_frequency);
	beta2 = getBeta(epsilon_r2,mu_r2,conductivity2,angular_frequency);
	//factor1 = Complex.Magnitude(wave_impedance2)*mu_r1/Complex.Magnitude(wave_impedance1)/mu_r2;
	factor1 = Math.sqrt(mu_r1*epsilon_r1/(mu_r2*epsilon_r2));
	    for(int i = 0; i<x.length; i++){
		theta1 = i * Math.PI/180.0;
		x[i] = theta1;
		//MARKER THETA2
		//theta2 = Math.asin(beta1*Math.sin(theta1)/beta2);
		theta2 = Math.asin(Math.sin(theta1)*factor1);
		//System.out.println("theta1 = "+theta1/Math.PI*180.0+"  theta2 = "+theta2*180.0/Math.PI+"   eta1 = "+wave_impedance1+"   eta2 = "+wave_impedance2);
		if(parallel){
			tmp = Complex.Divide(
			Complex.Multiply(wave_impedance2,2*Math.cos(theta1)),
			    Complex.Add(
				Complex.Multiply(wave_impedance2,Math.cos(theta2)),
				Complex.Multiply(wave_impedance1,Math.cos(theta1))
			    )
			);
	    }
		else {
		        tmp = Complex.Divide(
		        Complex.Multiply(wave_impedance2,2.0*Math.cos(theta1)),
			    Complex.Add(
				Complex.Multiply(wave_impedance2,Math.cos(theta1)),
				Complex.Multiply(wave_impedance1,Math.cos(theta2))
			    )
			);
		}
		tmp2 =  tmp.Magnitude()*tmp.Magnitude()
			*Complex.Magnitude(Complex.Divide(wave_impedance1,wave_impedance2))
			*Math.abs(Math.cos(theta2)/Math.cos(theta1));
		if((mu_r1*epsilon_r1) <= (mu_r2*epsilon_r2)){
		    array[i] = tmp2;
		}
		else{
		    //factor2 = Math.sin(theta1)*Math.sin(theta1)*factor1*factor1;
		    factor2 = Math.sin(theta1)*factor1;
		    if(factor2 > 1.0){
			array[i] = 0.0;
		    }
		    else{   
			array[i] = tmp2;
		    }
		}
	    }// end of loop in i
    }
}//End of EMF class


/** Bibliography
    [1] Pozar, David M.
	Microwave Engineering, 1990, Addison-Wesley Inc, ISBN 0-201-50418-9
    [2] Nannapaneni Narayana Rao
	Elements of Engineering Electromagnetics, Fourth Edition
	ISBN 0-13-948746-8, Prentice Hall Inc.
    [3] Liang Chi Shen / Jin Au Kong
	Applied Electromagnetism, Third Edition
	PWS Foundations in Engineering Series, ISBN 0-534-94722-0
 **/