far.c

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

/*****************************************************************************
 *  The following codes compute the far fields based on the surface current  *
 *  obtained by emap5 			                                         *
 *									                       *
 *  Authors:  Dr. Mohammad W. Ali,   Yun  Ji,  Dr. Todd. H. Hubing           *
 *  Version:  1.1							                       *
 *  Last updated:  Dec 16, 1998   				                       *
 *****************************************************************************/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "complex.c"
#include "alloc.c"
#include "util.c"


#ifndef M_PI  
#define M_PI 3.1415926
#endif

#define TotQuadPoint 7


/***********************  Function Prototype  *************************/
void   GetGaussQuadData();
void   ComputeGaussQuadPoint(int, int *, double *);
void   ReadInputFile();
void   ComputeNonSingul(int *, double, double, double, double, double);
void   ComputeFarEField(double *, double, double, double, double);
double Sign(int);


/***********************    GLOBAL VARIABLES   ************************/

int
   *InnerEdgeStat,
   *BoundEdgeStat,
   **PlusTrngleDetect,  *PlusTrngleIndex,
   **MinusTrngleDetect, *MinusTrngleIndex,
   **TrngleNode, **TrngleEdge, *BoundEdgeStat,
   **GlobalEdgeEnds,
   **TetEdge, **TetNode,
   TotEdgeNum,
   TotTetElement,
   TotNodeNum,
   TotTrngleNum,
   TotBoundEdgeNum,
   DielBoundEdgeNum,
   TotInnerEdgeNum,
   TotISourceEdgeNum,
   HybrdBoundEdgeNum;


double 
    WaveNumber, WaveLength, ImpeDance, FreeSpaceVel, AbsPermeable, AbsPermitt,
    OperateFreq,
    Wght[7], Qpnt[7][3], 
    RadDisObserv,
    **NodeCord;

complex
    InoPQ, IxiPQ, IetaPQ, IzetaPQ,  
    CFactor,
    CFactor1,  
    *Epsilon,
    *MVector,
    *JVector;

FILE *InF1,*InF2,*OutF;



main(int argc, char **argv)
{
    int IEdge, JEdge;
    double ThetaObserv1, ThetaObserv2, ThetaStep, 
           PhiObserv1, PhiObserv2, PhiStep, 
           ThetaObserv, PhiObserv,
           sth,sph,cth,cph, ObservPoint[3];  

    if (argc!=4)
    {
       fprintf(stderr,"Usage:far mesh SurfaceFieldFile OutputFile\n");
       exit(1);
    }

    if ((InF1  = fopen(argv[1],"r"))==NULL)
    {
       fprintf(stderr,"Can't open MeshFile:%s\n",argv[1]);
       exit(1);
    }

    if ((InF2 = fopen(argv[2],"r"))==NULL) 
    {
       fprintf(stderr,"Can't open SurfaceFieldFile:%s\n",argv[2]);
       exit(1) ;
    }

    if ((OutF = fopen(argv[3],"w"))==NULL)
    {
       fprintf(stderr,"Can't open OutputFile:%s\n",argv[3]);
       exit(1) ;
    }   

 
    GetGaussQuadData();
    ReadInputFile();
  
    FreeSpaceVel  = 3.0E+08;
    AbsPermeable  = 1.25663706144E-06;
    AbsPermitt    = 8.8542E-12;
    WaveLength = (FreeSpaceVel/(OperateFreq*1.0E+06));
    WaveNumber = 2.0*M_PI /WaveLength;
    ImpeDance  = WaveNumber/(2.0*M_PI*OperateFreq*1.0E+06*AbsPermitt);
    CFactor = Real_Mul(-1.0, COMplex_Cmplx(0.0,(WaveNumber*ImpeDance)));
           

    /* calculate the far field observed at (RadDisObserv,Theta,Phi)  */

    printf("Enter the observing distance(in wavelength, >=20): \t");
    scanf("%lf",&RadDisObserv);
    RadDisObserv *= WaveLength;

    printf("\nEnter Theta interval:\n\tStart Angle(in degree,>=0): \t");
    scanf("%lf", &ThetaObserv1);
    printf("\tEnd Angle(in degree, <=180): \t");
    scanf("%lf", &ThetaObserv2);
    printf("\tTheta step (!=0):\t");
    scanf("%lf", &ThetaStep);
 
    if (ThetaStep==0.0) 
    {
        fprintf(stderr, "Step can't be zero\n");
        exit(1);
    } 

    printf("\nEnter Phi interval:\n\tStart Angle(in degree, >=0): \t");
    scanf("%lf", &PhiObserv1);
    printf("\tEnd Angle(in degree,<=180): \t");
    scanf("%lf", &PhiObserv2);
    printf("\tPhi step(!=0):\t");
    scanf("%lf", &PhiStep);
 
    if (PhiStep==0.0) 
    {
        fprintf(stderr, "Step can't be zero\n");
        exit(1);
    }
     
    printf("\nPlease wait ......\n");
    fprintf(OutF, "Phi\t\tTheta\t\tE_Phi\t\tE_Theta\n");

    for(IEdge=0; IEdge<=(PhiObserv2-PhiObserv1)/PhiStep; IEdge++)
    {
        PhiObserv = PhiObserv1+ IEdge*PhiStep ;
        for(JEdge=0;JEdge<=(ThetaObserv2-ThetaObserv1)/ThetaStep;JEdge++)
        {
            ThetaObserv = ThetaObserv1+ JEdge*ThetaStep; 
            fprintf(OutF, "%lf\t%lf\t",  PhiObserv, ThetaObserv );  
            sth = sin(ThetaObserv * (M_PI/180.0));
            cth = cos(ThetaObserv * (M_PI/180.0));
            sph = sin(PhiObserv   * (M_PI/180.0));
            cph = cos(PhiObserv   * (M_PI/180.0));
            ObservPoint[0] = RadDisObserv * (sth * cph); 
            ObservPoint[1] = RadDisObserv * (sth * sph); 
            ObservPoint[2] = RadDisObserv * ( cth);
            ComputeFarEField(ObservPoint,cth,cph,sth,sph);
        }   /* for(JEdge */
    }  /* for(IEdge */

}  /* end of main()  */



/******************* ReadInputFile() *******************/
void ReadInputFile()
{
    int i,j ; 
    char buff[100], c;
 
    /* read the global information table */
    fscanf(InF1, " %c", &c);
    if (c!='#') 
    {
       fprintf(stderr, "Check your global information table\n");
       exit(1);
    }
    fgets(buff, sizeof(buff)-1, InF1);

    fscanf(InF1,"%d",&TotEdgeNum);
    fscanf(InF1,"%d",&TotTetElement);
    fscanf(InF1,"%d",&TotNodeNum);
    fscanf(InF1,"%d",&TotTrngleNum);
    fscanf(InF1,"%d",&TotBoundEdgeNum);
    fscanf(InF1,"%d",&TotInnerEdgeNum);
    fscanf(InF1,"%d",&DielBoundEdgeNum);
    fscanf(InF1,"%d",&HybrdBoundEdgeNum);

    /* read the global node table */
    fscanf(InF1, " %c", &c);
    if (c!='#') 
    {
       fprintf(stderr, "Check your global node table\n");
       exit(1);
    }
    fgets(buff, sizeof(buff)-1, InF1);
    NodeCord = double_Matrix(TotNodeNum,3);
    for(i=0;i<TotNodeNum;i++) 
        for(j=0;j<=2;j++)
            fscanf(InF1,"%lf",&NodeCord[i][j]);

    /* read the global edge table */
    fscanf(InF1, " %c", &c);
    if (c!='#') 
    {
       fprintf(stderr, "Check your global edge table\n");
       exit(1);
    }
    fgets(buff, sizeof(buff)-1, InF1);
    GlobalEdgeEnds= INT_Matrix(TotEdgeNum,6);
    for(i=0;i<TotEdgeNum;i++)
        for(j=0;j<6;j++)
            fscanf(InF1,"%d",&GlobalEdgeEnds[i][j]);


    /* read the triangle table */
    fscanf(InF1, " %c", &c);
    if ( c!='#' ) 
    {
       fprintf(stderr, "Check your triangle table\n");
       exit(1);
    }
    fgets(buff, sizeof(buff)-1, InF1);
    TrngleNode = INT_Matrix(TotTrngleNum,3);
    TrngleEdge = INT_Matrix(TotTrngleNum,3);
    for(i=0;i<TotTrngleNum;i++)
    {
        for(j=0;j<=2;j++) fscanf(InF1,"%d",&TrngleNode[i][j]);
        for(j=0;j<=2;j++) fscanf(InF1,"%d",&TrngleEdge[i][j]);
    }
 
  
    /* read tetrahedron table */
    fscanf(InF1, " %c", &c);
    if (c!='#') 
    {
       fprintf(stderr, "Check your tetrahedron table\n");
       exit(1);
    }
    fgets(buff, sizeof(buff)-1, InF1);
    TetNode = INT_Matrix(TotTetElement,4);
    Epsilon =CMPLX_Vector(TotTetElement);
    for(i=0;i<TotTetElement;i++)
    {
        for(j=0;j<=3;j++)
            fscanf(InF1,"%d",&TetNode[i][j]);
        fscanf(InF1,"%lf",&Epsilon[i].x);
        fscanf(InF1,"%lf",&Epsilon[i].y);
    }

 
    TetEdge = INT_Matrix(TotTetElement,6);
    for(i=0;i<TotTetElement;i++) for(j=0;j<=5;j++)
        fscanf(InF1,"%d",&TetEdge[i][j]);


    /* read plus/minus triangle table */    
    fscanf(InF1, " %c", &c);
    if (c!='#') 
    {
       fprintf(stderr, "Check your plus/minus triangle table\n");
       exit(1);
    }
    fgets(buff, sizeof(buff)-1, InF1);
    PlusTrngleDetect = INT_Matrix(TotEdgeNum,3);
    PlusTrngleIndex = INT_Vector(TotEdgeNum);
    MinusTrngleDetect = INT_Matrix(TotEdgeNum,3);
    MinusTrngleIndex = INT_Vector(TotEdgeNum);
    for(i=0;i<TotEdgeNum;i++)
    {
        fscanf(InF1,"%d",&PlusTrngleIndex[i]);
        for(j=0;j<PlusTrngleIndex[i];j++)
            fscanf(InF1,"%d",&PlusTrngleDetect[i][j]);
        fscanf(InF1,"%d",&MinusTrngleIndex[i]);