JSOC_Documentation: proj/mag/nlfff/apps/rebin.c Source File

00001 /*
00002  * Rebin magnetogram and field cube
00003  * Adapted from Wiegelmann's rebin.c & rebin2d.c
00004  * Adapted by X. Sun
00005  * Version 1.0 Mar 22 2010
00006  */
00007 
00008 /* Rebin magnetogram to 2 ^ (level - 1) smaller
00009  * Include level = 1 case
00010  * Adapted from rebin2d.c
00011  */
00012 void rebin2d(double *Bx, double *By, double *Bz, 
00013              double *Bx2, double *By2, double *Bz2, 
00014              int nx, int ny, int level)
00015 {
00016     int i, ix, iy;
00017     int i2, ix2, iy2, nx2, ny2;
00018     int k, kx, ky;
00019     int d;
00020     double dd;
00021     
00022     d = 1;
00023     for (k = level; k > 1; k--) { d *= 2; }
00024     dd = 1.0 * d * d;
00025     
00026     nx2 = nx / d; ny2 = ny / d;
00027  printf("rebin2d"); fflush(stdout);
00028     for (iy2 = 0; iy2 < ny2; iy2++)
00029     for (ix2 = 0; ix2 < nx2; ix2++)
00030     {
00031         i2 = ny2 * ix2 + iy2;
00032         ix = ix2 * d;
00033         iy = iy2 * d;
00034         i = ny * ix + iy;
00035         Bx2[i2] = 0.0; By2[i2] = 0.0; Bz2[i2] = 0.0;
00036         for (ky = 0; ky < d; ky++)
00037         for (kx = 0; kx < d; kx++)
00038         {
00039             Bx2[i2] += Bx[i + kx * ny + ky];
00040             By2[i2] += By[i + kx * ny + ky];
00041             Bz2[i2] += Bz[i + kx * ny + ky];
00042         }
00043         Bx2[i2] /= dd; By2[i2] /= dd; Bz2[i2] /= dd;
00044     }
00045      printf("rebin2d_done"); fflush(stdout);
00046 }
00047 
00048 /* Rebin 3d field cube to twice big in every dimension
00049  * Adapted from rebin.c
00050  */
00051 void rebin3d(double *Bx, double *By, double *Bz, 
00052              double *Bx2, double *By2, double *Bz2, 
00053              int nx, int ny, int nz)
00054 {
00055     int nynz = ny * nz, nxnynz = nx * ny * nz;
00056     int i, ix, iy, iz;
00057     int i2, ix2, iy2, iz2;
00058     int i3, ix3, iy3, iz3;
00059     
00060     for(ix2 = 0; ix2 < 2 * nx; ix2++)
00061     for(iy2 = 0; iy2 < 2 * ny; iy2++)
00062     for(iz2 = 0; iz2 < 2 * nz; iz2++)
00063     {
00064         ix = ix2 / 2;
00065         iy = iy2 / 2;
00066         iz = iz2 / 2;
00067         i = ix * nynz + iy * nz + iz;
00068         i2 = ix2 * 4 * nynz + iy2 * 2 * nz + iz2;
00069 
00070         ix3 = ix2 % 2;  /* For Trilinear interpolation */
00071         iy3 = iy2 % 2;
00072         iz3 = iz2 % 2;
00073         i3 = i;
00074         if (ix < nx - 1) { i3 = i3 + ix3 * nynz; }
00075         if (iy < ny - 1) { i3 = i3 + iy3 * nz; }
00076         if (iz < nz - 1) { i3 = i3 + iz3; }
00077         Bx2[i2] = 0.5 * (Bx[i] + Bx[i3]);
00078         By2[i2] = 0.5 * (By[i] + By[i3]);
00079         Bz2[i2] = 0.5 * (Bz[i] + Bz[i3]);
00080     }
00081 }