home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / SIMSEC.C < prev    next >
Encoding:
Text File  |  1984-07-27  |  4.9 KB  |  208 lines
  1.    simsec(a,b,x,y,tol,numrt)
  2.  
  3.       /*given two second-order simultaneous equations in two variables,*/
  4.       /*converts to fourth-order equation and solves for both variables*/
  5.       /*equations of the form                                          */
  6.       /*a(1)*x**2 + a(2)*x*y + a(3)*y**2 + a(4)*x + a(5)*y + a(6) = 0  */
  7.       /*b(1)*x**2 + b(2)*x*y + b(3)*y**2 + b(4)*x + b(5)*y + b(6) = 0  */
  8.  
  9.       int *numrt;
  10.       float a[],b[],x[],y[],tol;
  11.  
  12.     {
  13.       int i,iroot,j,jj,k,kk,l,n,nortx,norty,nroot;
  14.       float c[5],d[5],w1[6],w2[6],rt[12],root[12],zero;
  15.       float dum1[2],dum2[2],dum3[2],dum4[2],c3sqd,c4sqd;
  16.       extern double fabs();
  17.  
  18.  
  19.       zero = tol/10.;
  20.       if((fabs(a[0]) - zero) > 0.)
  21.       {
  22.       if((fabs(b[0]) - zero) >0) goto nq;
  23.       for( i = 0; i <= 5; i++)
  24.        {
  25.        w1[i] = b[i];
  26.        w2[i] = a[i];
  27.        }
  28.       }
  29.        else
  30.        for(i = 0; i <= 5; i++)
  31.        {
  32.         w1[i] = a[i];
  33.         w2[i] = b[i];
  34.        }
  35.  
  36.       if((fabs(w1[1]) - zero) <= 0. && (fabs(w1[3]) - zero) <= 0.)
  37.       {
  38.        d[0] = w1[2];
  39.        d[1] = w1[4];
  40.        d[2] = w1[5];
  41.        quadrt(d,dum1,dum2,tol,&norty);
  42.        rt[0] = dum1[0];
  43.        rt[1] = dum1[1];
  44.        rt[2] = dum2[0];
  45.        rt[3] = dum2[1];
  46.        goto l25;
  47.       }
  48.       c[0] = w1[2];
  49.       c[1] = w1[4];
  50.       c[2] = w1[5];
  51.       c[3] = -w1[1];
  52.       c[4] = -w1[3];
  53.       goto l12;
  54.  
  55. nq:   for(i = 0; i <= 5; i++)
  56.       {
  57.        w1[i] = a[i];
  58.        w2[i] = b[i];
  59.       }
  60.  
  61.       for(i = 0; i <= 1; i++)
  62.       {
  63.        j = (2*(i+1));
  64.        c[i] = w2[0]*w1[j]-w1[0]*w2[j];
  65.        j --;
  66.        c[i+3] = w1[0]*w2[j]-w2[0]*w1[j];
  67.       }
  68.       c[2] = w2[0]*w1[5] - w1[0]*w2[5];
  69. l12:  c3sqd = c[3]*c[3];
  70.       c4sqd = c[4]*c[4];
  71.       d[0] = w2[0]*c[0]*c[0] + w2[1]*c[0]*c[3] + w2[2]*c3sqd;
  72.  
  73.       d[1] = 2.*w2[0]*c[0]*c[1] + w2[1]*(c[0]*c[4] + c[1]*c[3]) + 2.*
  74.         w2[2]*c[3]*c[4] + w2[3]*c[0]*c[3] + w2[4]*c3sqd;
  75.  
  76.       d[2] = w2[0]*(2.*c[0]*c[2]+c[1]*c[1])+w2[1]*(c[1]*c[4]+c[2]*
  77.         c[3]) + w2[2]*c4sqd + w2[3]*(c[0]*c[4] + c[1]*c[3]) + 2.
  78.         *w2[4]*c[3]*c[4] + w2[5]*c3sqd;
  79.  
  80.       d[3] = 2.*w2[0]*c[1]*c[2] + w2[1]*c[2]*c[4] + w2[3]*
  81.         (c[1]*c[4] + c[2]*c[3]) + w2[4]*c4sqd + 2.*w2[5]*c[3]*c[4];
  82.  
  83.       d[4] = w2[0]*c[2]*c[2] + w2[3]*c[2]*c[4] + w2[5]*c4sqd;
  84.  
  85.       biquad(d,dum1,dum2,dum3,dum4,tol,&nroot);
  86.       rt[0] = dum1[0];
  87.       rt[1] = dum1[1];
  88.       rt[2] = dum2[0];
  89.       rt[3] = dum2[1];
  90.       rt[4] = dum3[0];
  91.       rt[5] = dum3[1];
  92.       rt[6] = dum4[0];
  93.       rt[7] = dum4[1];
  94.       *numrt = nroot;
  95.  
  96. l25:  if((*numrt - 1) < 0) return;
  97.       k = 0;
  98.       n = *numrt;
  99.       for(i = 1; i <= n; i++)
  100.       {
  101.        kk = (2 * i) -1;
  102.        if((fabs(rt[kk]) - zero) <= 0.)
  103.        {
  104.         rt[k] = rt[kk - 1];
  105.         k++;
  106.        }
  107.       }
  108.  
  109.       if((k - 1) < 0)
  110.       {
  111.        *numrt = 0;
  112.        return;
  113.       }
  114.       if((k - 1) > 0)
  115.       {
  116.        norty = 0;
  117.        kk = k - 2;
  118.        for(i = 0; i <= kk; i++)
  119.        {
  120.         for(j = i+1; j <= k-1; j++)
  121.         {
  122.          if((fabs(rt[i]-rt[j])-zero) <= 0.) goto l15;
  123.         }
  124.         norty = norty + 1;
  125.         rt[norty-1] = rt[i];
  126. l15:   ;
  127.        }
  128.        norty = norty + 1;
  129.        rt[norty-1] = rt[k-1];
  130.       }
  131.        else norty = 1;
  132.       nortx = 0;
  133.  
  134.       for(j = 0; j <= norty-1; j++)
  135.       {
  136.        k = 0;
  137.        for(i = 0; i <= 6; i += 6)
  138.        {
  139.         if(i == 0)
  140.         {
  141.          c[0] = w1[0];
  142.          c[1] = w1[1] *rt[j] + w1[3];
  143.          c[2] = w1[2] * rt[j]*rt[j] + w1[4]*rt[j]+w1[5];
  144.         }
  145.          else
  146.          {
  147.           c[0] = w2[0];
  148.           c[1] = w2[1] *rt[j] + w2[3];
  149.           c[2] = w2[2] * rt[j]*rt[j] + w2[4]*rt[j]+w2[5];
  150.          }
  151.         quadrt(c,dum1,dum2,tol,&iroot);
  152.         root[k]   = dum1[0];
  153.         root[k+1] = dum1[1];
  154.         root[k+2] = dum2[0];
  155.         root[k+3] = dum2[1];
  156.  
  157.         if(i == 0)
  158.         {
  159.         *numrt = iroot;
  160.         k = k + iroot*2;
  161.         }
  162.        }
  163.        if(*numrt <= 0 )
  164.        {
  165.        rt[j] = 99999.;
  166.        goto l27;
  167.        }
  168.        if(iroot  <= 0 )
  169.        {
  170.         rt[j] = 99999.;
  171.         goto l27;
  172.        }
  173.       jj = *numrt * 2;
  174.       l = jj + 2;
  175.       kk = (iroot + *numrt) * 2;
  176.  
  177.       for(i = 1; i <= jj-1; i+=2)
  178.       {
  179.        for(k = l-1; k <= kk-1; k+=2)
  180.        {
  181.         if((fabs(root[k]) - zero) <= 0.)
  182.  
  183.         {
  184.          if((fabs(root[i]) - zero) > 0.) goto l30;
  185.          if((fabs(root[k-1] - root[i-1]) - zero) <= 0.) goto l34;
  186.         }
  187.        }
  188.        goto l30;
  189. l34:   nortx = nortx + 1;
  190.        x[nortx - 1] = root[k - 1];
  191.        goto l27;
  192. l30:  ;
  193.       }
  194.       rt[j] = 99999.;
  195. l27:  ;
  196.       }
  197.       *numrt = 0;
  198.       for(i = 0; i <= norty-1; i++)
  199.       {
  200.        if((rt[i] - 99999.) != 0.)
  201.        {
  202.         *numrt = *numrt + 1;
  203.         y[*numrt-1] = rt[i];
  204.        }
  205.       }
  206.       return;
  207.     }
  208.