home *** CD-ROM | disk | FTP | other *** search
/ Los Alamos National Laboratory / LANL_CD.ISO / software / compres / src / us_conv2.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_lp_rbc(x1, x2, h, npt, m)
  8. int npt, m;
  9. float *x1, *h;
  10. register float *x2;
  11. {
  12.     int nptr, nptl, flag, nptl2, nptr2;
  13.     register int m2, m21;
  14.     float *max1, *max2, *max3, *hmax2, *h1, *x1p1, *x1m1;
  15.     register float *xptr, *hmax1, *hptr, *ptr1, *ptr2;
  16.  
  17.     nptl = m2 = m/2;
  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.     x1p1 = x1 + 1;
  25.     x1m1 = x1 - 1;
  26.  
  27.     flag = ( m2 % 2 ) ? 1 : 0;
  28.     max1 = x2 + nptl - 1;
  29.     max2 = x2 + 2*npt - nptr - 1;
  30.     max3 = ( flag ) ? x2 + 2*npt - 1 : x2 + 2*npt;
  31.  
  32. /*________________________________________________________________*/
  33.  
  34. /*  Type 2 reflection  */
  35.  
  36.     ptr1 = x1 + nptl2;
  37.     ptr2 = ( flag ) ? x1m1 + nptr2 : x1 + nptr2;
  38.     if(!flag) {
  39.         nptl2--;
  40.         nptr2++;
  41.     }
  42.     hptr = ( flag ) ? h1 : h;
  43.     hmax1 = h + nptl--;
  44.     *x2 = 0.;
  45.     while(hptr < hmax1) {
  46.         *x2 += *hptr * (*ptr1-- + *ptr2--);
  47.         hptr += 2;
  48.     }
  49.     *x2++ += *hptr * *ptr2;
  50.  
  51.  
  52.     if(flag) {
  53.         ptr1 = x1 + nptl2--;
  54.         ptr2 = x1 + nptr2++;
  55.         hptr = h;
  56.         hmax1 = h1 + nptl--;
  57.         *x2 = 0.;
  58.         while(hptr < hmax1) {
  59.             *x2 += *hptr * (*ptr1-- + *ptr2--);
  60.             hptr += 2;
  61.     }
  62.         ptr1 = x1p1;
  63.         while(hptr < hmax2) {
  64.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  65.             hptr += 2;
  66.     }
  67.         if( !flag )
  68.             *x2++ += *hptr * *ptr2;
  69.         else
  70.             x2++;
  71.       }
  72.  
  73.  
  74.     while(x2 < max1) {
  75.         ptr1 = x1 + nptl2;
  76.         ptr2 = x1m1 + nptr2;
  77.         hptr = h1;
  78.         hmax1 = (hptr = h1) + nptl--;
  79.         *x2 = 0.;
  80.         while(hptr < hmax1) {
  81.             *x2 += *hptr * (*ptr1-- + *ptr2--);
  82.             hptr += 2;
  83.     }
  84.  
  85.         ptr1 = x1p1;
  86.         while(hptr < hmax2) {
  87.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  88.             hptr += 2;
  89.     }
  90.         if( flag )
  91.             *x2++ += *hptr * *ptr2;
  92.         else
  93.             x2++;
  94.  
  95.  
  96.         ptr1 = x1 + nptl2--;
  97.         ptr2 = x1 + nptr2++;
  98.         hptr = h;
  99.         hmax1 = h1 + nptl--;
  100.         *x2 = 0.;
  101.         while(hptr < hmax1) {
  102.             *x2 += *hptr * (*ptr1-- + *ptr2--);
  103.             hptr += 2;
  104.     }
  105.         ptr1 = x1p1;
  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.         hptr = h1;
  122.         hmax1 = hptr + m21;
  123.         *x2 = 0.;
  124.         while(hptr < hmax1) {
  125.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  126.             hptr += 2;
  127.     }
  128.         if( flag )
  129.             *x2++ += *hptr * *ptr1;
  130.         else
  131.             x2++;
  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.     ptr2 = (ptr1 = xptr) + m21;
  147.     hmax1 = (hptr  = h1) + m21;
  148.     *x2 = 0.;
  149.     while(hptr < hmax1) {
  150.         *x2 += *hptr * (*ptr1++ + *ptr2--);
  151.         hptr += 2;
  152.     }
  153.     if( flag )
  154.         *x2++ += *hptr * *ptr1;
  155.     else
  156.         x2++;
  157.  
  158. /*________________________________________________________________*/
  159.  
  160. /*  Type 1 reflection  */
  161.  
  162.     nptr2 = (nptr = 1) / 2;
  163.     nptl2 = (m - nptr) /2;
  164.     while(x2 < max3) {
  165.         hmax1 = (hptr = h) + nptr++;
  166.         ptr1 = x1 + npt - nptl2--;
  167.         ptr2 = x1 + npt - ++nptr2;
  168.         *x2 = 0.;
  169.         while(hptr < hmax1) {
  170.             *x2 += *hptr * (*ptr1++ + *ptr2++);
  171.             hptr += 2;
  172.     }
  173.  
  174.         ptr2 = x1m1 + npt;
  175.         while(hptr < hmax2) {
  176.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  177.             hptr += 2;
  178.     }
  179.         if( !flag )
  180.             *x2++ += *hptr * *ptr2;
  181.         else
  182.             x2++;
  183.  
  184.  
  185.         hptr = h1;
  186.         hmax1 = h + nptr++;
  187.         ptr1 = x1 + npt - nptl2;
  188.         ptr2 = x1 + npt - nptr2;
  189.         *x2 = 0.;
  190.         while(hptr < hmax1) {
  191.             *x2 += *hptr * (*ptr1++ + *ptr2++);
  192.             hptr += 2;
  193.     }
  194.         ptr2 = x1m1 + npt;
  195.         while(hptr < hmax2) {
  196.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  197.             hptr += 2;
  198.     }
  199.  
  200.         if( flag )
  201.             *x2++ += *hptr * *ptr2;
  202.         else
  203.             x2++;
  204.     }
  205.  
  206.  
  207.     if( flag ) {
  208.         hmax1 = (hptr = h) + nptr;
  209.         ptr1 = x1 + npt - nptl2--;
  210.         ptr2 = x1 + npt - ++nptr2;
  211.         *x2 = 0.;
  212.         while(hptr < hmax1) {
  213.             *x2 += *hptr * (*ptr1++ + *ptr2++);
  214.             hptr += 2;
  215.     }
  216.  
  217.         ptr2 = x1;
  218.         while(hptr < hmax2) {
  219.             *x2 += *hptr * (*ptr1++ + *ptr2--);
  220.             hptr += 2;
  221.     }
  222.         if( !flag )
  223.             *x2 += *hptr * *ptr2;
  224.     }
  225. }
  226.  
  227.