src/mesh/geom.cc

   1 /*****************************************************************************/
   2 /*              8888888         88888888         88888888                    */
   3 /*            8                8                8                            */
   4 /*           8                 8                8                            */
   5 /*           8                  88888888         88888888                    */
   6 /*           8      8888                8                8                   */
   7 /*            8       8                 8                8                   */
   8 /*              888888         888888888        888888888                    */
   9 /*                                                                           */
  10 /*       A Two-Dimensional General Purpose Semiconductor Simulator.          */
  11 /*                                                                           */
  12 /*  GSS 0.4x                                                                 */
  13 /*  Last update: May 29, 2006                                                */
  14 /*                                                                           */
  15 /*  Gong Ding                                                                */
  16 /*  gdiso@ustc.edu                                                           */
  17 /*  NINT, No.69 P.O.Box, Xi'an City, China                                   */
  18 /*                                                                           */
  19 /*****************************************************************************/
  20
  21 #include <math.h>
  22 #include <stdio.h>
  23 #include "geom.h"
  24
  25
  26 /* ----------------------------------------------------------------------------
  27  * CompAngle:  This function computes the angle between three points given
  28  * by coordinates (x1,y1), (x2,y2) and (x3,y3).  The angle is measured from
  29  * vector 12 to 32.clock order
  30  */
  31 double CompAngle(double x1, double y1,double x2, double y2,double x3, double y3)
  32 {
  33   double xl,yl;                         /* coordinates of vector 12 */
  34   double xr,yr;                         /* coordinates of vector 32 */
  35   double r;                             /* temporary number */
  36   double angle;                         /* angle between vectors 12 & 32 */
  37   double PI=3.14159265358979323846;
  38   xl = x1 - x2;  yl = y1 - y2;
  39   xr = x3 - x2;  yr = y3 - y2;
  40
  41   r = (xl*xr+yl*yr)/(sqrt(xl*xl+yl*yl)*sqrt(xr*xr+yr*yr));
  42   if (r > 1.0)
  43     r =  1.0;
  44   else if (r < -1.0) r = -1.0;
  45   angle = acos(r);
  46   if (xl*yr-xr*yl > 0)
  47     angle = 2*PI-angle;
  48   return(angle);
  49
  50 } /* CompAngle */
  51
  52 /* ----------------------------------------------------------------------------
  53  * CircleCenter:  This function determines the center of the circle that
  54  * passes through points (x1,y1), (x2,y2) and (x3,y3).  If the points are
  55  * coplaner then 1 is returned otherwise 0 is
  56  * returned.
  57  */
  58 int CircleCenter(double x1, double y1, double x2, double y2,double x3, double y3,
  59                   double &xc, double &yc)
  60 {
  61   double xl,yl;                         /* coordinates of vector 12 */
  62   double xr,yr;                         /* coordinates of vector 32 */
  63   double rl,rr;                         /* ??? */
  64   double det;                           /* determinent of coordinate matrix */
  65
  66   xl = x1 - x2;  yl = y1 - y2;
  67   xr = x3 - x2;  yr = y3 - y2;
  68
  69   det = xr * yl - xl * yr;
  70   if (det == 0.0)
  71   {
  72     xc = yc = 0.0;
  73     return 1;
  74   }
  75
  76   rl = (xl * (x1 + x2) + yl * (y1 + y2)) / det;
  77   rr = (xr * (x3 + x2) + yr * (y3 + y2)) / det;
  78
  79   xc = -0.5 * (rl * yr - rr * yl);
  80   yc = -0.5 * (rr * xl - rl * xr);
  81   return 0;
  82 } /* CircleCenter */
  83
  84
  85 /* ----------------------------------------------------------------------------
  86  * MirrorPoint:  This function computes the MirrorPoint(xg,yg) of (x,y) with the
  87  * line determined by (x1,y1) and (x2,y2).
  88  */
  89 void MirrorPoint(double x,double y,double x1,double y1,double x2,double y2,double &xg,double &yg)
  90 {
  91   double r,ang1,ang2,angle;
  92
  93
  94   double r1 = sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y));
  95   double r2 = sqrt((x2-x)*(x2-x)+(y2-y)*(y2-y));
  96   if(r1>r2)
  97   {
  98     r = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
  99     ang1 = CompAngle(x2,y2,x1,y1,x,y);
 100     ang2 = CompAngle(x2,y2,x1,y1,x1+r,y1);
 101     angle = ang1 + ang2;
 102     xg = x1 + r1*cos(angle);
 103     yg = y1 + r1*sin(angle);
 104   }
 105   else
 106   {
 107     r = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 108     ang1 = CompAngle(x1,y1,x2,y2,x,y);
 109     ang2 = CompAngle(x1,y1,x2,y2,x2+r,y2);
 110     angle = ang1 + ang2;
 111     xg = x2 + r2*cos(angle);
 112     yg = y2 + r2*sin(angle);
 113
 114   }
 115 }/* CompAngle */
 116
 117
 118 /* ----------------------------------------------------------------------------
 119  * DistancePointLine:  This function computes the distance between point (x,y) and
 120  * line determined by (x1,y1) and (x2,y2).
 121  */
 122 double DistancePointLine(double x1, double y1,double x2, double y2,double x, double y)
 123 {
 124   double LineMag = Distance(x1,y1,x2,y2);
 125   double U = (((x-x1)*(x2-x1))+((y-y1)*(y2-y1))) / (LineMag*LineMag);
 126
 127   double IntersectionX = x1 + U*(x2-x1);
 128   double IntersectionY = y1 + U*(y2-y1);
 129
 130   return  Distance( x, y, IntersectionX, IntersectionY);
 131 }/* DistancePointLine */
 132
 133
 134 /* ----------------------------------------------------------------------------
 135  * Distance:  This function computes the distance between (x1,y1) and (x2,y2)
 136  * line determined by (x1,y1) and (x2,y2).
 137  */
 138 double Distance(double x1, double y1,double x2,double y2)
 139 {
 140   if(x1==x2 && y1==y2) return 0;
 141   return sqrt((x1-x2)*(x1-x2) +(y1-y2)*(y1-y2));
 142 }/* Distance */
 143
 144
 145 /* ----------------------------------------------------------------------------
 146  * TriArea:  This function computes the area of a triangle determined by (x1,y1)
 147  * (x2,y2) and (x3,y3).
 148  */
 149 double TriArea(double x1, double y1,double x2, double y2,double x3, double y3)
 150 {
 151   double d1,d2,d3,d,r,xc,yc;
 152   d1=Distance(x1,y1,x2,y2);
 153   d2=Distance(x1,y1,x3,y3);
 154   d3=Distance(x2,y2,x3,y3);
 155   d=0.5*(d1+d2+d3);
 156   return sqrt(d*(d-d1)*(d-d2)*(d-d3));
 157 }/* TriArea */
 158
 159
 160 /* ----------------------------------------------------------------------------
 161  * TriAngle:  This function compute the max angle of a triangle
 162  *
 163  */
 164 double MaxTriAngle(double x1, double y1, double x2, double y2,double x3, double y3)
 165 {
 166   double a,b,c;
 167   a=Distance(x1,y1,x2,y2);
 168   b=Distance(x1,y1,x3,y3);
 169   c=Distance(x2,y2,x3,y3);
 170   if(a>b&&a>c) return acos((b*b+c*c-a*a)/(2*b*c));
 171   if(b>a&&b>c) return acos((a*a+c*c-b*b)/(2*a*c));
 172   if(c>b&&c>a) return acos((a*a+b*b-c*c)/(2*a*b));
 173   return 0; //prevent warning of complier
 174 } /* TriAngle */
 175
 176
 177 /* ----------------------------------------------------------------------------
 178  * RaySegmentIntersectTest:  This function do 2D Ray-Segment Intersection test
 179  *
 180  */
 181 IntersectResult RaySegmentIntersectTest(double begin_x, double begin_y, double ex, double ey, double x1,double y1,double x2,double y2)
 182 {
 183   double end_x =  begin_x + ex/sqrt(ex*ex+ey*ey+1e-20);
 184   double end_y =  begin_y + ey/sqrt(ex*ex+ey*ey+1e-20);
 185   double denom  = ((y2 - y1)*(end_x - begin_x)) -((x2 - x1)*(end_y - begin_y));
 186   double nume_a = ((x2 - x1)*(begin_y - y1)) - ((y2 - y1)*(begin_x - x1));
 187   double nume_b = ((end_x - begin_x)*(begin_y - y1)) - ((end_y - begin_y)*(begin_x - x1));
 188   if(denom == 0.0)
 189   {
 190     if(nume_a == 0.0 && nume_b == 0.0)
 191     {
 192       return COINCIDENT;
 193     }
 194     return PARALLEL;
 195   }
 196
 197   double ua = nume_a / denom;
 198   double ub = nume_b / denom;
 199
 200   if(ub >= 0.0 && ub <= 1.0)
 201   {
 202     // Get the intersection point.
 203     double intersection_x = begin_x + ua*(end_x - begin_x);
 204     double intersection_y = begin_y + ua*(end_y - begin_y);
 205     return INTERESECTING;
 206   }
 207
 208   return NOT_INTERESECTING;
 209 }
 210