complex.c

三维矢量有限元-矩量法电磁场分析程序。 EMAP5 is a full-wave electromagnetic field solver that combines the method of m

/***************************************************************************
 *  The following functions support complex variable operations 	         *
 *									                     *
 *  Authors:  Dr. Mohammad W. Ali,   Yun  Ji,  Dr. Todd. H. Hubing         *
 *  Version:  1.1							                     *
 *  Last updated:  Dec 16, 1998   				                     *
 ***************************************************************************/



/********************   Define "complex" Type   **************************/

typedef struct { 
    double x;
    double y; 
} complex;


/******************   Function  Prototype Defintion  ***********************/

double  Aimag(complex);
double  COMplex_Abs(complex);
complex COMplex_Add(complex, complex);
complex COMplex_Add2(complex, complex, complex);
complex COMplex_Cmplx(double, double);
complex COMplex_Conjg(complex);
complex COMplex_Div(complex, complex);
complex COMplex_Expon(double, complex);
complex COMplex_Mul(complex, complex);
complex COMplex_Null();
complex COMplex_Sub(complex, complex);
void    Null_CMPLXMatrix(complex **, int, int);
void    Null_CMPLXVector(complex *,int);
double  Real(complex );
complex Real_Div(complex, double);
complex Real_Mul(double, complex);


/***************************************************************************
Prototype:    double    Aimag( complex z) 
Description:     To get the imaginary part of a complex number 
Input value: 
    complex z  --- a complex number 
Return value:    the imaginary part of z.   
****************************************************************************/
double Aimag( complex z)
{
    return z.y;
}



/***************************************************************************
Prototype:    double    COMplex_Abs( complex z) 
Description:     To get the Euclid normal of a complex number 
Input value: 
    complex z  --- a complex number 
Return value:    Assume z=x+j y, return sqrt(x^2+y^2).   
****************************************************************************/
double COMplex_Abs( complex z)
{
 
    /* "Numerical recipe" P171 */
  
    if(fabs(z.x)<1e-30 && fabs(z.y)<1e-30)
           return 1e-30;
    else if(fabs(z.x)>fabs(z.y))  
             return (fabs(z.x)*sqrt(1.0+pow(z.y/z.x, 2.0)) );
    else return (fabs(z.y)*sqrt(1.0+pow(z.x/z.y, 2.0)) );
    
}

 
    
/***************************************************************************
Prototype:    complex    COMplex_Add( complex a, complex b) 
Description:     To add two complex numbers. 
Input value: 
    complex a ,b ---  two complex number 
Return value:     the complex sum of a and b   
****************************************************************************/
complex COMplex_Add(complex a, complex b)
{
    return COMplex_Cmplx(a.x + b.x,a.y + b.y);
}




/***************************************************************************
Prototype:    complex    COMplex_Add2(complex a, complex b, complex c) 
Description:     To add three complex variables 
Input value: 
    complex a, b, c --- three complex numbers 
Return value:     the complex sum of a, b and c  
****************************************************************************/
complex COMplex_Add2(complex a, complex b, complex c)
{
    return COMplex_Cmplx(a.x + b.x + c.x,a.y + b.y + c.y);
}




/****************************************************************************
Prototype:     complex    COMplex_Cmplx(double x , double y); 
Description:     To initialize  a complex number. 
Input value: 
    double x -- the real part 
    double y -- the image part 
Return value: a complex value equal to x+jy, where j is the imaginary unit. 
****************************************************************************/
complex COMplex_Cmplx(double x, double y)
{ 
    complex z; 

    z.x = x; 
    z.y = y; 

    return z;
}




/***************************************************************************
Prototype:    complex    COMplex_Conjg( complex z) 
Description:     To get the conjungate of a complex number. 
Input value: 
    complex z ---- a complex number 
Return value:     The conjugate of z.  
****************************************************************************/
complex COMplex_Conjg(complex z)
{ 
    return COMplex_Cmplx(z.x, -z.y);
}



/***************************************************************************
Prototype:    complex    Real_Div(complex z, double a) 
Description:    To divide a complex number by a real number. 
Input value: 
    complex z --- a complex number 
    double a --- a double number 
Return value:     the complex number of z / a   
****************************************************************************/
complex COMplex_Div(complex a, complex b)
{ 
    double f1,f2;
    
    if(fabs((double)b.x)>fabs((double)b.y)) {
        f1=b.y/b.x;
        f2=b.x+b.y*f1;
        return COMplex_Cmplx( (a.x+a.y*f1)/f2, (a.y-a.x*f1)/f2);
    }
       else {
            f1=b.x/b.y;
            f2=b.y+b.x*f1;
            return COMplex_Cmplx( (a.y+a.x*f1)/f2, (-a.x+a.y*f1)/f2);
       }  
    
}      



/***************************************************************************
Prototype:    complex    COMplex_Expon(double a, complex z) 
Description:     To compute the exponential  function z**a. 
Input value: 
    complex z  --- the base number 
    double a  --- the exponential number 
Return value:     the exponential value of z**a. 
Subroutines called:     COMplex_Cmplx()   
****************************************************************************/
complex COMplex_Expon(double a, complex z)
{ 
    double R = exp((double)(a*z.x));
    return COMplex_Cmplx(R*cos(a*z.y),R*sin(a*z.y));
}



/***************************************************************************
Prototype:    complex    COMplex_Mul(complex a, complex b) 
Description:     To multiple two complex variables 
Input value: 
    complex a, b  --- two complex numbers 
Return value:     the complex product of a and b. 
Subroutines called:     COMplex_Cmplx()   
****************************************************************************/
complex COMplex_Mul(complex a, complex b)
{ 
    return COMplex_Cmplx(a.x*b.x - a.y*b.y,a.x*b.y + a.y*b.x);
}





/***************************************************************************
Prototype:    complex    COMplex_Null() 
Description:     To  initialize a complex to zero 
Input value:    none 
Return value:     a complex number equals to zero. 
****************************************************************************/
complex COMplex_Null()
{ 
    complex z; 
    
    z.x=0.0; 
    z.y=0.0; 
  
    return z; 
}
 



/***************************************************************************
Prototype:    complex    COMplex_Sub(complex a, complex b) 
Description:     To subtract complex b from complex a 
Input value: 
    complex a, b --- two complex numbers 
Return value:    the complex number of  a-b   
****************************************************************************/
complex COMplex_Sub(complex a, complex b)
{   
    return COMplex_Cmplx(a.x - b.x,a.y - b.y);
}




/***************************************************************************
Prototype:    void     Null_CMPLXVector( complex *Vec, int LenX) 
Description    To initialize a vector of complex type to zero 
Input value: 
    complex *Vec  --- pointer to a complex vector 
    int LenX --- the length of the vector 
Return value:    none   
****************************************************************************/
void Null_CMPLXVector(complex *Vec, int LenX)
{
    int ii;
    
    for(ii=0;ii<LenX;ii++)
        Vec[ii] = COMplex_Null();

}




/***************************************************************************
Prototype:    double    Null_CMPLXMAtrix(complex **Vec, int LenX, int LenY) 
Description    To initialize a matrix of complex type to zero. 
Input value: 
    complex **Vec  --- pointer to a matrix of complex type 
    int LenX  ---  row length 
    int LenY  --- column length 
Return value:     none   
****************************************************************************/
void Null_CMPLXMatrix(complex **Vec, int LenX, int LenY)
{
    int ii,jj;

    for(ii=0;ii<LenX;ii++)
        for(jj=0;jj<LenY;jj++)
            Vec[ii][jj] = COMplex_Null();

}



/***************************************************************************
Prototype:    double    Real(complex z) 
Description    To get the real part of a complex number 
Input value: 
    complex z  --- a complex number 
Return value:     the real part of z.   
****************************************************************************/
double Real(complex z)
{  
    return z.x;
}



/***************************************************************************
Prototype:    complex    Real_Div(complex z, double a) 
Description:    To divide a complex number by a real number. 
Input value: 
    complex z --- a complex number 
    double a --- a double number 
Return value:     the complex number of z / a   
****************************************************************************/
complex Real_Div(complex z, double a)
{
    return COMplex_Cmplx(z.x/a,z.y/a);
}
  


/***************************************************************************
Prototype:    complex    Real_Mul(double a, complex z) 
Description:     To get the product of a real number and a complex number 
Input value: 
    double a --- a real number 
    complex z  --- a complex number 
Return value:     complex product of a and z 
Subroutines called:     COMplex_Cmplx() 
****************************************************************************/
complex Real_Mul(double a, complex z)
{
   return COMplex_Cmplx(a*z.x,a*z.y);
}