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

00001 /*
00002  * Magnetogram preprocessing
00003  * Adapted from Wiegelmann's prepro.c
00004  * Takes in bz0 magnetogram
00005  * Computes and outputs processed mag
00006  * Adapted by X. Sun
00007  * Version 1.0 Mar 21 2010
00008  *  Parameters needs to be tweaked
00009  */
00010 
00011 
00012 
00013 // mu1 = 1.0; mu2 = 1.0; mu3 = 0.001; mu4 = 0.01; mu5 = 0.0;
00014 // maxit = 10000;
00015 
00016 void preproc(double *Bx, double *By, double *Bz, 
00017              int nx, int ny, 
00018              double mu1, double mu2, double mu3, double mu4, double mu5, 
00019              int maxit, int verb)
00020 {
00021     double *x, *y, *r;
00022 //  double *Bx, *By, *Bz;
00023 //  double *Bxold, *Byold, *Bzold;
00024     double *Bxo, *Byo, *Bzo;
00025     double *Bx1, *By1, *Bz1;
00026 //  double *Hx, *Hy, *w;
00027     double *LapBx, *LapBy, *LapBz;
00028     double *term3a, *term3b, *term3c;
00029     double *term4a, *term4b, *term4c;
00030 //  char leer[5];
00031     double Bave, dt;    //, dumm1, fac1, dt0
00032 //  double mu1, mu2, mu3, mu4, mu5;
00033     double dx, dy, dxdy;
00034     double Emag, ihelp; //, Emag2, ihelp2
00035     double dL, L, L1, L2, L12, L3, L4, L5;
00036     double oldL12, oldL3, oldL4, oldL5, oldL;
00037     double dL12, dL3, dL4, dL5;
00038     double term1a, term1b, term1c, term1c1, term1c2; 
00039     double term2a, term2b, term2c;
00040     double term2a1, term2a2, term2b1, term2b2, term2c1, term2c2;
00041     double term5;
00042 //  double Lvec[5];
00043     int nxny;   //, nx, ny; 
00044     int i, ix, iy, it;
00045     int i1, i2, ix1, iy1;
00046     int i3, i4, ix2, iy2;
00047 
00048     nxny = nx * ny;
00049 
00050 //  Bx =(double *) calloc(nxny, sizeof(double));
00051 //  By =(double *) calloc(nxny, sizeof(double));
00052 //  Bz =(double *) calloc(nxny, sizeof(double));
00053     LapBx = (double *) calloc(nxny, sizeof(double));
00054     LapBy = (double *) calloc(nxny, sizeof(double));
00055     LapBz = (double *) calloc(nxny, sizeof(double));
00056     Bx1 = (double *) calloc(nxny, sizeof(double));
00057     By1 = (double *) calloc(nxny, sizeof(double));
00058     Bz1 = (double *) calloc(nxny, sizeof(double));
00059     Bxo = (double *) calloc(nxny, sizeof(double));
00060     Byo = (double *) calloc(nxny, sizeof(double));
00061     Bzo = (double *) calloc(nxny, sizeof(double));
00062 //  Bxold = (double *) calloc(nxny, sizeof(double));
00063 //  Byold = (double *) calloc(nxny, sizeof(double));
00064 //  Bzold = (double *) calloc(nxny, sizeof(double));
00065     x = (double *) calloc(nx, sizeof(double));
00066     y = (double *) calloc(ny, sizeof(double));
00067     r = (double *) calloc(nxny, sizeof(double));
00068     term3a = (double *) calloc(nxny, sizeof(double));
00069     term3b = (double *) calloc(nxny, sizeof(double));
00070     term3c = (double *) calloc(nxny, sizeof(double));
00071     term4a = (double *) calloc(nxny, sizeof(double));
00072     term4b = (double *) calloc(nxny, sizeof(double));
00073     term4c = (double *) calloc(nxny, sizeof(double));
00074 /*
00075     if(mu5 !=0.0) {
00076         printf("\n Using Halpha images for preprocessing");
00077         Hx =(DOUBLE *) calloc(nxny, sizeof(DOUBLE));
00078         Hy =(DOUBLE *) calloc(nxny, sizeof(DOUBLE));
00079         w =(DOUBLE *) calloc(nxny, sizeof(DOUBLE));
00080         if((streamw = fopen("halpha.dat","r")) == NULL) { printf("\n Error halpha.dat "); exit(1);}
00081         for(iy=0;iy<ny;iy++)
00082         for(ix=0;ix<nx;ix++)
00083         {
00084         i=ny*ix+iy; 
00085         fscanf(streamw,"%lf", &dummy1);
00086         Hx[i]=dummy1;
00087         fscanf(streamw,"%lf", &dummy1);
00088         Hy[i]=dummy1;
00089         }
00090         fclose(streamw);
00091         printf("\n Halpha image (Hx,Hy) halpha.dat loaded ");
00092         if((streamw = fopen("w.dat","r")) == NULL) { printf("\n Error w.dat "); exit(1);}
00093         for(iy=0;iy<ny;iy++)
00094         for(ix=0;ix<nx;ix++)
00095         {
00096         i=ny*ix+iy; 
00097         fscanf(streamw,"%lf", &dummy1);
00098         w[i]=dummy1;
00099         }
00100         fclose(streamw);
00101         printf("\n weighting function w.dat loaded ");
00102     }
00103  */
00104 
00105     if (verb) printf("Raw vectormagnetogram loaded \n");
00106     Bave = 0.0;
00107 #ifdef _OPENMP
00108 #pragma omp parallel for private (i,ix,iy) reduction(+:Bave)
00109 #endif
00110     for (ix = 0; ix < nx; ix++)
00111     for (iy = 0; iy < ny; iy++) {
00112         i = ny * ix + iy;
00113         Bave = Bave + sqrt(Bx[i] * Bx[i] + By[i] * By[i] + Bz[i] * Bz[i]);
00114     }
00115     Bave = Bave / nx / ny;
00116 //  printf("\n Bave= %lf",Bave);
00117 
00118 #ifdef _OPENMP
00119 #pragma omp parallel for private (i)
00120 #endif
00121     for (i = 0; i < nxny; i++) {
00122         Bx[i] = Bx1[i] = Bxo[i] = Bx[i] / Bave;
00123         By[i] = By1[i] = Byo[i] = By[i] / Bave;
00124         Bz[i] = Bz1[i] = Bzo[i] = Bz[i] / Bave;
00125     }
00126     dx = 1.0 / (nx - 1); dy = 1.0 / (ny - 1);
00127     dxdy = dx * dy;
00128     for (ix = 0; ix < nx; ix++) x[ix] = ix * dx;
00129     for (iy = 0; iy < ny; iy++) y[iy] = iy * dy;
00130 
00131 #ifdef _OPENMP
00132 #pragma omp parallel for private (i,ix,iy)
00133 #endif
00134     for (ix = 0; ix < nx; ix++)
00135     for (iy = 0; iy < ny; iy++) {
00136         i = ny * ix + iy;
00137         r[i] = sqrt(x[ix] * x[ix] + y[iy] * y[iy]);
00138     }
00139 
00140     /* Do the Preprocessing */
00141 //  fac1 = 10.0;
00142     dt = 0.01;
00143 //  printf("\n dt changed to:");
00144 //  printf("\n dt= %lf", dt);
00145     Emag = ihelp = 0.0;
00146 #ifdef _OPENMP
00147 #pragma omp parallel for private (i) reduction(+:Emag) reduction(+:ihelp)
00148 #endif
00149     for (i = 0; i < nxny; i++) {
00150         Emag = Emag + (Bx[i] * Bx[i] + By[i] * By[i] + Bz[i] * Bz[i]);
00151         ihelp = ihelp + (r[i] * (Bx[i] * Bx[i] + By[i] * By[i] + Bz[i] * Bz[i]));
00152     }
00153 //  printf("\n Emag= %lf, ihelp= %lf",Emag,ihelp);
00154 
00155     /*************************************/
00156     /********** Start Iteration **********/
00157     /*************************************/
00158     L = 10.0; dL=1.0; oldL=10.0;
00159     it = -1;
00160 //  printf("\n L12, L3, L4 are multiplied with 1.0e6");
00161     while(it < maxit && dL > 0.0001 && dt > 1.0e-6)
00162     {
00163         it++;
00164         term1a = term1b = term1c = term1c1 = term1c2 = 0.0;
00165         term2a = term2b = term2c = 0.0;
00166         term2a1 = term2b1 = term2c1 = 0.0;
00167         term2a2 = term2b2 = term2c2 = 0.0;
00168         L1 = L2 = L3 = L4 = L5 = 0.0;
00169 #ifdef _OPENMP
00170 #pragma omp parallel for private (i,ix,iy) reduction(+:term1a) reduction(+:term1b) reduction(+:term1c1) reduction(+:term1c2) reduction(+:term2a1) reduction(+:term2a2) reduction(+:term2b1) reduction(+:term2b2) reduction(+:term2c1) reduction(+:term2c2) reduction(+:L3)
00171 #endif
00172         for (ix = 0; ix < nx; ix++)
00173         for (iy = 0; iy < ny; iy++)
00174         {
00175             i = ny * ix + iy;
00176             /********** Term 1 Force **********/
00177             term1a = term1a + Bx[i] * Bz[i];
00178             term1b = term1b + By[i] * Bz[i];
00179             term1c1 = term1c1 + Bz[i] * Bz[i];
00180             term1c2 = term1c2 + Bx[i] * Bx[i] + By[i] * By[i];
00181             /********** Term 2 Torque **********/
00182             term2a1 = term2a1 + x[ix] * Bz[i] * Bz[i];
00183             term2a2 = term2a2 + x[ix] * (Bx[i] * Bx[i] + By[i] * By[i]);
00184             term2b1 = term2b1 + y[iy] * Bz[i] * Bz[i];
00185             term2b2 = term2b2 + y[iy] * (Bx[i] * Bx[i] + By[i] * By[i]);
00186             term2c1 = term2c1 + y[iy] * Bx[i] * Bz[i];
00187             term2c2 = term2c2 + x[ix] * By[i] * Bz[i];
00188             /********** Term 3 Data **********/
00189             term3a[i] = (Bx[i] - Bxo[i]) * (Bx[i] - Bxo[i]) / Emag;
00190             term3b[i] = (By[i] - Byo[i]) * (By[i] - Byo[i]) / Emag;
00191             term3c[i] = (Bz[i] - Bzo[i]) * (Bz[i] - Bzo[i]) / Emag;
00192             L3 = L3 + term3a[i] + term3b[i] + term3c[i];
00193         }
00194         term1a = term1a / Emag;
00195         term1b = term1b / Emag;
00196         term1c = (term1c1 - term1c2) / Emag;
00197         term2a = (term2a1 - term2a2) / ihelp;
00198         term2b = (term2b1 - term2b2) / ihelp;
00199         term2c = (term2c1 - term2c2) / ihelp;
00200         L1 = term1a * term1a + term1b * term1b + term1c * term1c;
00201         L2 = term2a * term2a + term2b * term2b + term2c * term2c;
00202         L12 = L1 + L2;
00203         L3 = L3 / Emag;
00204         /********** Laplace Bx,By,Bz **********/
00205 #ifdef _OPENMP
00206 #pragma omp parallel for private (i,ix,iy,ix1,ix2,iy1,iy2,i1,i2,i3,i4) reduction(+:L4)
00207 #endif
00208         for (ix = 0; ix < nx; ix++)
00209         for (iy = 0; iy < ny; iy++)
00210         {
00211             i = ny * ix + iy;
00212             ix1 = ix + 1; if (ix1 == nx) ix1 = 0;
00213             ix2 = ix - 1; if (ix2 == -1) ix2 = nx - 1;
00214             iy1 = iy + 1; if (iy1 == ny) iy1 = 0;
00215             iy2 = iy - 1; if (iy2 == -1) iy2 = ny - 1;
00216             i1 = ny * ix1 + iy;
00217             i2 = ny * ix2 + iy;
00218             i3 = ny * ix + iy1;
00219             i4 = ny * ix + iy2;
00220             LapBx[i] = -4 * Bx[i] + Bx[i1] + Bx[i2] + Bx[i3] + Bx[i4];
00221             LapBy[i] = -4 * By[i] + By[i1] + By[i2] + By[i3] + By[i4];
00222             LapBz[i] = -4 * Bz[i] + Bz[i1] + Bz[i2] + Bz[i3] + Bz[i4];
00223             L4 = L4 + LapBx[i] * LapBx[i] + LapBy[i] * LapBy[i] + LapBz[i] * LapBz[i];
00224         }
00225         L4 = dxdy * L4 / Emag;
00226 //        L = L1 + L2 + L3 + L4 + L5;//??
00227         L = mu1 * L1 + mu2 * L2 + mu3 * L3 + mu4 * L4 + mu5 * L5;
00228         /********** Term 4 Smoothing **********/
00229 #ifdef _OPENMP
00230 #pragma omp parallel for private (i,ix,iy,ix1,ix2,iy1,iy2,i1,i2,i3,i4)
00231 #endif
00232         for(ix = 0; ix < nx; ix++)
00233         for(iy = 0; iy < ny; iy++)
00234         {
00235             i = ny * ix + iy;
00236             ix1 = ix + 1; if (ix1 == nx) ix1 = 0;
00237             ix2 = ix - 1; if (ix2 == -1) ix2 = nx - 1;
00238             iy1 = iy + 1; if (iy1 == ny) iy1 = 0;
00239             iy2 = iy - 1; if (iy2 == -1) iy2 = ny - 1;
00240             i1 = ny * ix1 + iy;
00241             i2 = ny * ix2 + iy;
00242             i3 = ny * ix + iy1;
00243             i4 = ny * ix + iy2;
00244             term4a[i] = - 8 * LapBx[i] + 2 * LapBx[i1] 
00245                         + 2 * LapBx[i2] + 2 * LapBx[i3] + 2 * LapBx[i4];
00246             term4b[i] = - 8 * LapBy[i] + 2 * LapBy[i1] 
00247                         + 2 * LapBy[i2] + 2 * LapBy[i3] + 2 * LapBy[i4];
00248             term4c[i] = - 8 * LapBz[i] + 2 * LapBz[i1] 
00249                         + 2 * LapBz[i2] + 2 * LapBz[i3] + 2 * LapBz[i4];
00250         }
00251         if (it > 0) {
00252             dL12 = dL3 = dL4 = dL5 = 0.0;
00253             if (L12 != 0.0) dL12 = fabs(L12 - oldL12) / L12;
00254             if (L3 != 0.0) dL3 = fabs(L3 - oldL3) / L3;
00255             if (L4 != 0.0) dL4 = fabs(L4 - oldL4) / L4;
00256             if (L5 != 0.0) dL5 = fabs(L5 - oldL5) / L5;
00257             /* dL=dL12+dL3+dL4; */
00258             if (L != 0.0) dL = (oldL - L) / L;
00259         }
00260         /*if(it%50==0)*/
00261 //        if ((it % 100 == 0 || it > 3355) && verb == 2) {
00262         if (it % 100 == 0 && verb == 2) {
00263             printf("%i L12= %lf, L3= %lf, L4= %lf, L=%lf, dL= %lf\n", 
00264             it, 1.0e6*L12, 1.0e6*L3, 1.0e6*L4, 1.0e6*L, dL);
00265         }
00266         if(oldL < L && it > 0) {
00267             dt = dt / 2;
00268             if (verb == 2) printf("dt reduced: %i, dt=%lf, dL=%lf, oldL=%lf, L=%lf\n",
00269             it, dt, dL, 1.0e6*oldL, 1.0e6*L);
00270             dL = 1.0;
00271             oldL = L;
00272         } else {
00273             if (it > 0) dt = dt * 1.001;
00274             oldL12 = L12; oldL3 = L3;
00275             oldL4 = L4; oldL5 = L5;
00276             oldL = L;
00277         }
00278         /********** Iterate for new B **********/
00279 #ifdef _OPENMP
00280 #pragma omp parallel for private (i,ix,iy)
00281 #endif
00282         for (ix = 0; ix < nx; ix++)
00283         for (iy = 0; iy < ny; iy++)
00284         {
00285             i = ny * ix + iy;
00286             Bx1[i] = Bx[i] + dt * 
00287                      (mu1 * (-2 * term1a * Bz[i] + 4 * term1c * Bx[i])
00288                       - mu3 * 2 * (Bx[i] - Bxo[i])
00289                       - mu4 * term4a[i]
00290                       + mu2 * 
00291                           (4 * term2a * x[ix] * Bx[i] + 4 * term2b * y[iy] * Bx[i] 
00292                            - 2 * term2c * y[iy] * Bz[i]));
00293             By1[i] = By[i] + dt * 
00294                      (mu1 * (-2 * term1b * Bz[i] + 4 * term1c * By[i])
00295                       - mu3 * 2 * (By[i] - Byo[i])
00296                       - mu4 * term4b[i]
00297                       + mu2 * 
00298                           (4 * term2a * x[ix] * By[i] + 4 * term2b * y[iy] * By[i]
00299                            + 2 * term2c * x[ix] * Bz[i]));
00300             Bz1[i] = Bz[i] + dt * (- mu3 * 2 * (Bz[i] - Bzo[i]) - mu4 * term4c[i]);
00301         }
00302         /********** Halpha terms ******/
00303 /*
00304         if(mu5 !=0.0)
00305         {
00306             term5=L5=0.0;
00307             #ifdef _OPENMP
00308             #pragma omp parallel for private (i,ix,iy) reduction(+:L5)
00309             #endif
00310             for(ix=0;ix<nx;ix++)
00311             for(iy=0;iy<ny;iy++)
00312             {
00313                 i=ny*ix+iy;
00314                 term5=Bx[i]*Hy[i]-By[i]*Hx[i];
00315                 Bx1[i]=Bx1[i]-dt*w[i]*mu5*term5*Hy[i];
00316                 By1[i]=By1[i]+dt*w[i]*mu5*term5*Hx[i];
00317                 L5=L5+term5*term5;
00318             }
00319             L5=dxdy*L5/Emag;
00320             if(it%50==0) {printf(" L5= %lf",1.0e6*L5);}
00321         }
00322  */
00323         /********** Update B **********/
00324 #ifdef _OPENMP
00325 #pragma omp parallel for private (i)
00326 #endif
00327         for (i = 0; i < nxny; i++) {
00328             Bx[i] = Bx1[i];
00329             By[i] = By1[i];
00330             Bz[i] = Bz1[i];
00331         }
00332 
00333     } /*End of Iteration*/
00334     
00335     if (verb) printf("End of Iteration\n");
00336     if (verb == 2) {
00337         printf("mu3= %lf \t mu4= %lf \t mu5= %lf\n", mu3, mu4, mu5);
00338         printf("%i L12= %lf, L3= %lf, L4= %lf, L=%lf, dL= %lf\n", 
00339         it, 1.0e6*L12, 1.0e6*L3, 1.0e6*L4, 1.0e6*L,dL);
00340     }
00341     
00342     /* Output data */
00343 #ifdef _OPENMP
00344 #pragma omp parallel for private (i)
00345 #endif
00346     for (i = 0; i < nxny; i++) {
00347         Bx[i] = Bx[i] * Bave;
00348         By[i] = By[i] * Bave;
00349         Bz[i] = Bz[i] * Bave;
00350     }
00351     
00352     /* Clean up */
00353     free(LapBx); free(LapBy); free(LapBz);
00354     free(Bx1); free(By1); free(Bz1);
00355     free(Bxo); free(Byo); free(Bzo);
00356 //  free(Bxold); free(Byold); free(Bzold);
00357     free(x); free(y); free(r);
00358     free(term3a); free(term3b); free(term3c);
00359     free(term4a); free(term4b); free(term4c);
00360 
00361 }