home *** CD-ROM | disk | FTP | other *** search
/ Los Alamos National Laboratory / LANL_CD.ISO / software / compres / src / us_conv1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-02-11  |  4.9 KB  |  227 lines
  1. /****************************************************************
  2.  
  3. COPYRIGHT (C) 1992 UNIVERSITY OF CALIFORNIA
  4.  
  5. ***************************************************************/
  6.  
  7. us_conv_hp_rbc(x1, x2, h, npt, m)
  8. int npt, m;
  9. register float *x1, *h;
  10. register float *x2;
  11. {
  12.     int nptr, nptl, flag, nptl2, nptr2, npt1;
  13.     register int m2, m21;
  14.     float *max1, *max2, *max3, *hmax2, *h1, *x1m1;
  15.     register float *xptr, *hmax1, *hptr, *ptr1, *ptr2;
  16.  
  17.     nptl = (m2 = m/2) + 1;
  18.     nptr = m - nptl;
  19.     nptl2 = nptl/2;
  20.     nptr2 = nptr/2;
  21.     m21 = m2 - 1;
  22.     h1 = h + 1;
  23.     hmax2 = h + m2;
  24.     x1m1 = x1 - 1;
  25.  
  26.     flag = ( m2 % 2 ) ? 1 : 0;
  27.     max1 = x2 + nptl - 1;
  28.     max2 = x2 + 2*npt - nptr - 1;
  29.     max3 = ( !flag ) ? x2 + 2*npt - 1 : x2 + 2*npt;
  30.  
  31. /*________________________________________________________________*/
  32.  
  33. /*  Type 1 reflection  */
  34.  
  35.     if(!flag) {
  36.         nptl--;
  37.         hptr = h1;
  38.         hmax1 = h + nptr++;
  39.         ptr1 = x1m1 + nptl2;
  40.         ptr2 = x1m1 + nptr2;
  41.         *x2 = 0.;
  42.         while(hptr < hmax1) {
  43.             *x2 += *hptr * (*ptr1-- + *ptr2--);
  44.             hptr += 2;
  45.         }
  46.         ptr2 = x1;
  47.         while(hptr < hmax2) {
  48.             *x2 += *hptr * (*ptr1-- + *ptr2++);
  49.             hptr += 2;
  50.       }
  51.         x2++;
  52.     }
  53.  
  54.  
  55.     nptl--;
  56.     hmax1 = (hptr = h) + nptr++ - 1;
  57.     ptr1 = x1m1 + nptl2--;
  58.     ptr2 = x1 + nptr2++;
  59.     *x2 = 0.;
  60.     while(hptr < hmax1) {
  61.         *x2 += *hptr * (*ptr1-- + *ptr2--);
  62.         hptr += 2;
  63.     }
  64.  
  65.     ptr2 = x1;
  66.     while(hptr < hmax2) {
  67.         *x2 += *hptr * (*ptr1-- + *ptr2++);
  68.         hptr += 2;
  69.     }
  70.     if( !flag )
  71.         *x2++ += *hptr * *ptr2;
  72.     else
  73.         x2++;
  74.  
  75.  
  76.     while(x2 < max1) {
  77.         ptr1 = x1m1 + nptl2;
  78.         ptr2 = x1m1 + nptr2;
  79.         hptr = h1;
  80.         hmax1 = h + nptl--;
  81.         *x2 = 0.;
  82.         while(hptr < hmax1) {
  83.             *x2 += *hptr * (*ptr1-- + *ptr2--);
  84.             hptr += 2;
  85.     }
  86.         ptr1 = x1;
  87.         while(hptr < hmax2) {
  88.             *x2 += *hptr * (*ptr1-- + *ptr2--);
  89.             hptr += 2;
  90.     }
  91.         if( flag )
  92.             *x2++ += *hptr * *ptr2;
  93.         else
  94.             x2++;
  95.  
  96.  
  97.         hmax1 = (hptr = h) + nptl--;
  98.         ptr1 = x1m1 + nptl2--;
  99.         ptr2 = x1 + nptr2++;
  100.         *x2 = 0.;
  101.         while(hptr < hmax1) {
  102.             *x2 += *hptr * (*ptr1-- + *ptr2--);
  103.             hptr += 2;
  104.     }
  105.         ptr1 = x1;
  106.         while(hptr < hmax2) {
  107.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  108.             hptr += 2;
  109.     }
  110.         if( !flag )
  111.             *x2++ += *hptr * *ptr2;
  112.         else
  113.             x2++;
  114.     }
  115.  
  116. /*________________________________________________________________*/
  117.  
  118.     xptr = x1;
  119.     while(x2 < max2) {
  120.         ptr2 = (ptr1 = xptr) + m21;
  121.         hmax1 = (hptr = h1) + m21;
  122.         *x2 = 0.;
  123.         while(hptr < hmax1) {
  124.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  125.             hptr += 2;
  126.     }
  127.         if( flag )
  128.             *x2++ += *hptr * *ptr1;
  129.         else
  130.             x2++;
  131.  
  132.  
  133.         ptr2 = (ptr1 = xptr) + m2;
  134.         hmax1 = (hptr = h) + m2;
  135.         *x2 = 0.;
  136.         while(hptr < hmax1) {
  137.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  138.             hptr += 2;
  139.     }
  140.         if( !flag )
  141.             *x2++ += *hptr * *ptr1;
  142.         else
  143.             x2++;
  144.         xptr++;
  145.     }
  146.  
  147.  
  148.     ptr2 = (ptr1 = xptr) + m21;
  149.     hmax1 = (hptr = h1) + m21;
  150.     *x2 = 0.;
  151.     while(hptr < hmax1) {
  152.         *x2 += *hptr * (*ptr1++ + *ptr2--);
  153.         hptr += 2;
  154.     }
  155.     if( flag )
  156.         *x2++ += *hptr * *ptr1;
  157.     else
  158.         x2++;
  159.  
  160. /*________________________________________________________________*/
  161.  
  162. /*  Type 2 reflection  */
  163.  
  164.     npt1 = npt - 1;
  165.     nptr2 = (nptr = 1) / 2;
  166.     nptl2 = (m - nptr) /2;
  167.     while(x2 < max3) {
  168.         hmax1 = (hptr = h) + nptr++;
  169.         ptr1 = x1 + npt - nptl2--;
  170.         ptr2 = x1m1 + npt - ++nptr2;
  171.         *x2 = 0.;
  172.         while(hptr < hmax1) {
  173.             *x2 += *hptr * (*ptr1++ + *ptr2++);
  174.             hptr += 2;
  175.     }
  176.  
  177.         ptr2 = x1m1 + npt;
  178.         while(hptr < hmax2) {
  179.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  180.             hptr += 2;
  181.     }
  182.         if( !flag )
  183.             *x2++ += *hptr * *ptr2;
  184.         else
  185.             x2++;
  186.  
  187.  
  188.         hptr = h1;
  189.         hmax1 = h + nptr++;
  190.         ptr1 = x1 + npt - nptl2;
  191.         ptr2 = x1m1 + npt - nptr2;
  192.         *x2 = 0.;
  193.         while(hptr < hmax1) {
  194.             *x2 += *hptr * (*ptr1++ + *ptr2++);
  195.             hptr += 2;
  196.     }
  197.         ptr2 = x1m1 + npt;
  198.         while(hptr < hmax2) {
  199.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  200.             hptr += 2;
  201.     }
  202.         if( flag )
  203.             *x2++ += *hptr * *ptr2;
  204.         else
  205.             x2++;
  206.     }
  207.  
  208.     if( !flag ) {
  209.         nptl = m - nptr;
  210.         hmax1 = (hptr = h) + nptr;
  211.         ptr1 = x1 + npt - nptl/2;
  212.         ptr2 = x1m1 + npt - ++nptr2;
  213.         *x2 = 0.;
  214.         while(hptr < hmax1) {
  215.             *x2 += *hptr * (*ptr1++ + *ptr2++);
  216.             hptr += 2;
  217.     }
  218.  
  219.         ptr2 = x1 + npt1;
  220.         while(hptr < hmax2) {
  221.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  222.             hptr += 2;
  223.     }
  224.         *x2 += *hptr * *ptr2;
  225.     }
  226. }
  227.