00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include <math.h>
00035 #include <stdio.h>
00036 #include <stdlib.h>
00037 #include "../defs.h"
00038 #include "../athena.h"
00039 #include "../globals.h"
00040 #include "prototypes.h"
00041 #include "../prototypes.h"
00042
00043 #ifdef THIRD_ORDER_CHAR
00044 #ifdef SPECIAL_RELATIVITY
00045 #error : PPM reconstruction (order=3) cannot be used for special relativity.
00046 #endif
00047
00048 static Real **pW=NULL, **dWm=NULL, **Wim1h=NULL;
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 void lr_states(const GridS* pG, const Prim1DS W[], const Real Bxc[],
00068 const Real dt, const Real dx, const int il, const int iu,
00069 Prim1DS Wl[], Prim1DS Wr[], const int dir)
00070 {
00071 int i,n,m;
00072 Real lim_slope1,lim_slope2,qa,qb,qc,qx;
00073 Real ev [NWAVE],rem [NWAVE][NWAVE],lem [NWAVE][NWAVE];
00074 Real ev_ip1[NWAVE],rem_ip1[NWAVE][NWAVE],lem_ip1[NWAVE][NWAVE];
00075 Real dWc[NWAVE+NSCALARS],dWl[NWAVE+NSCALARS];
00076 Real dWr[NWAVE+NSCALARS],dWg[NWAVE+NSCALARS];
00077 Real dac[NWAVE+NSCALARS],dal[NWAVE+NSCALARS];
00078 Real dar[NWAVE+NSCALARS],dag[NWAVE+NSCALARS],da[NWAVE+NSCALARS];
00079 Real Wlv[NWAVE+NSCALARS],Wrv[NWAVE+NSCALARS];
00080 Real dW[NWAVE+NSCALARS],W6[NWAVE+NSCALARS];
00081 Real *pWl, *pWr;
00082 Real qx1,qx2;
00083
00084
00085 Real ql,qr,qxx1,qxx2,zc,zr,zl,q1,q2,gamma_curv;
00086 int hllallwave_flag = 0;
00087 const Real dtodx = dt/dx;
00088 #ifdef CYLINDRICAL
00089 const Real *r=pG->r, *ri=pG->ri;
00090 #endif
00091
00092
00093 for (n=0; n<NWAVE; n++) {
00094 for (m=0; m<NWAVE; m++) {
00095 rem[n][m] = 0.0;
00096 lem[n][m] = 0.0;
00097 rem_ip1[n][m] = 0.0;
00098 lem_ip1[n][m] = 0.0;
00099 }
00100 }
00101 for (i=il-3; i<=iu+3; i++) pW[i] = (Real*)&(W[i]);
00102
00103
00104
00105
00106
00107
00108 for (i=il-2; i<=il-1; i++) {
00109
00110 #ifdef HYDRO
00111 #ifdef ISOTHERMAL
00112 esys_prim_iso_hyd(W[i].d,W[i].Vx, ev,rem,lem);
00113 #else
00114 esys_prim_adb_hyd(W[i].d,W[i].Vx,Gamma*W[i].P,ev,rem,lem);
00115 #endif
00116 #endif
00117
00118 #ifdef MHD
00119 #ifdef ISOTHERMAL
00120 esys_prim_iso_mhd(
00121 W[i].d,W[i].Vx, Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00122 #else
00123 esys_prim_adb_mhd(
00124 W[i].d,W[i].Vx,Gamma*W[i].P,Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00125 #endif
00126 #endif
00127
00128
00129
00130
00131
00132 for (n=0; n<(NWAVE+NSCALARS); n++) {
00133 #ifdef CYLINDRICAL
00134 if (dir==1) {
00135 dWc[n] = (pW[i+1][n]*r[i+1] - pW[i-1][n]*r[i-1])/r[i];
00136 dWl[n] = (pW[i ][n]*r[i ] - pW[i-1][n]*r[i-1])/ri[i];
00137 dWr[n] = (pW[i+1][n]*r[i+1] - pW[i ][n]*r[i ])/ri[i+1];
00138 }
00139 else
00140 #endif
00141 {
00142 dWc[n] = pW[i+1][n] - pW[i-1][n];
00143 dWl[n] = pW[i ][n] - pW[i-1][n];
00144 dWr[n] = pW[i+1][n] - pW[i ][n];
00145 }
00146 if (dWl[n]*dWr[n] > 0.0) {
00147 dWg[n] = 2.0*dWl[n]*dWr[n]/(dWl[n]+dWr[n]);
00148 } else {
00149 dWg[n] = 0.0;
00150 }
00151 }
00152
00153
00154
00155
00156 for (n=0; n<NWAVE; n++) {
00157 dac[n] = lem[n][0]*dWc[0];
00158 dal[n] = lem[n][0]*dWl[0];
00159 dar[n] = lem[n][0]*dWr[0];
00160 dag[n] = lem[n][0]*dWg[0];
00161 for (m=1; m<NWAVE; m++) {
00162 dac[n] += lem[n][m]*dWc[m];
00163 dal[n] += lem[n][m]*dWl[m];
00164 dar[n] += lem[n][m]*dWr[m];
00165 dag[n] += lem[n][m]*dWg[m];
00166 }
00167 }
00168
00169
00170
00171
00172 #if (NSCALARS > 0)
00173 for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00174 dac[n] = dWc[n];
00175 dal[n] = dWl[n];
00176 dar[n] = dWr[n];
00177 dag[n] = dWg[n];
00178 }
00179 #endif
00180
00181
00182
00183
00184 for (n=0; n<(NWAVE+NSCALARS); n++) {
00185 da[n] = 0.0;
00186 if (dal[n]*dar[n] > 0.0) {
00187 lim_slope1 = MIN( fabs(dal[n]),fabs(dar[n]));
00188 lim_slope2 = MIN(0.5*fabs(dac[n]),fabs(dag[n]));
00189 da[n] = SIGN(dac[n])*MIN(2.0*lim_slope1,lim_slope2);
00190 }
00191 }
00192
00193
00194
00195
00196 for (n=0; n<NWAVE; n++) {
00197 dWm[i][n] = da[0]*rem[n][0];
00198 for (m=1; m<NWAVE; m++) {
00199 dWm[i][n] += da[m]*rem[n][m];
00200 }
00201 }
00202
00203 #if (NSCALARS > 0)
00204 for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00205 dWm[i][n] = da[n];
00206 }
00207 #endif
00208
00209
00210
00211
00212
00213
00214
00215
00216 #ifdef H_CORRECTION
00217
00218
00219
00220
00221
00222
00223
00224
00225 #endif
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 }
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260 for (n=0; n<(NWAVE+NSCALARS); n++) {
00261 #ifdef CYLINDRICAL
00262 if (dir==1) {
00263 Wim1h[il-1][n] = (0.5*(pW[il-1][n]*r[il-1] + pW[il-2][n]*r[il-2])
00264 - (dWm[il-1][n]*r[il-1] - dWm[il-2][n]*r[il-2])/6.0)/ri[il-1];
00265 }
00266 else
00267 #endif
00268 {
00269 Wim1h[il-1][n] =.5*(pW[il-1][n]+pW[il-2][n])-(dWm[il-1][n]-dWm[il-2][n])/6.;
00270 }
00271 }
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 for (i=il-1; i<=iu+1; i++) {
00286
00287
00288
00289
00290 #ifdef HYDRO
00291 #ifdef ISOTHERMAL
00292 esys_prim_iso_hyd(W[i+1].d,W[i+1].Vx, ev_ip1,rem_ip1,lem_ip1);
00293 #else
00294 esys_prim_adb_hyd(W[i+1].d,W[i+1].Vx,Gamma*W[i+1].P,ev_ip1,rem_ip1,lem_ip1);
00295 #endif
00296 #endif
00297
00298 #ifdef MHD
00299 #ifdef ISOTHERMAL
00300 esys_prim_iso_mhd(W[i+1].d,W[i+1].Vx,
00301 Bxc[i+1],W[i+1].By,W[i+1].Bz,ev_ip1,rem_ip1,lem_ip1);
00302 #else
00303 esys_prim_adb_mhd(W[i+1].d,W[i+1].Vx,Gamma*W[i+1].P,
00304 Bxc[i+1],W[i+1].By,W[i+1].Bz,ev_ip1,rem_ip1,lem_ip1);
00305 #endif
00306 #endif
00307
00308
00309
00310
00311 for (n=0; n<(NWAVE+NSCALARS); n++) {
00312 #ifdef CYLINDRICAL
00313 if (dir==1) {
00314 dWc[n] = (pW[i+2][n]*r[i+2] - pW[i ][n]*r[i ])/r[i+1];
00315 dWl[n] = (pW[i+1][n]*r[i+1] - pW[i ][n]*r[i ])/ri[i+1];
00316 dWr[n] = (pW[i+2][n]*r[i+2] - pW[i+1][n]*r[i+1])/ri[i+2];
00317 }
00318 else
00319 #endif
00320 {
00321 dWc[n] = pW[i+2][n] - pW[i ][n];
00322 dWl[n] = pW[i+1][n] - pW[i ][n];
00323 dWr[n] = pW[i+2][n] - pW[i+1][n];
00324 }
00325 if (dWl[n]*dWr[n] > 0.0) {
00326 dWg[n] = 2.0*dWl[n]*dWr[n]/(dWl[n]+dWr[n]);
00327 } else {
00328 dWg[n] = 0.0;
00329 }
00330 }
00331
00332
00333
00334
00335 for (n=0; n<NWAVE; n++) {
00336 dac[n] = lem_ip1[n][0]*dWc[0];
00337 dal[n] = lem_ip1[n][0]*dWl[0];
00338 dar[n] = lem_ip1[n][0]*dWr[0];
00339 dag[n] = lem_ip1[n][0]*dWg[0];
00340 for (m=1; m<NWAVE; m++) {
00341 dac[n] += lem_ip1[n][m]*dWc[m];
00342 dal[n] += lem_ip1[n][m]*dWl[m];
00343 dar[n] += lem_ip1[n][m]*dWr[m];
00344 dag[n] += lem_ip1[n][m]*dWg[m];
00345 }
00346 }
00347
00348
00349
00350
00351 #if (NSCALARS > 0)
00352 for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00353 dac[n] = dWc[n];
00354 dal[n] = dWl[n];
00355 dar[n] = dWr[n];
00356 dag[n] = dWg[n];
00357 }
00358 #endif
00359
00360
00361
00362
00363 for (n=0; n<(NWAVE+NSCALARS); n++) {
00364 da[n] = 0.0;
00365 if (dal[n]*dar[n] > 0.0) {
00366 lim_slope1 = MIN( fabs(dal[n]),fabs(dar[n]));
00367 lim_slope2 = MIN(0.5*fabs(dac[n]),fabs(dag[n]));
00368 da[n] = SIGN(dac[n])*MIN(2.0*lim_slope1,lim_slope2);
00369 }
00370 }
00371
00372
00373
00374
00375 for (n=0; n<NWAVE; n++) {
00376 dWm[i+1][n] = da[0]*rem_ip1[n][0];
00377 for (m=1; m<NWAVE; m++) {
00378 dWm[i+1][n] += da[m]*rem_ip1[n][m];
00379 }
00380 }
00381
00382 #if (NSCALARS > 0)
00383 for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00384 dWm[i+1][n] = da[n];
00385 }
00386 #endif
00387
00388
00389
00390
00391
00392
00393
00394
00395 #ifdef H_CORRECTION
00396
00397
00398
00399
00400
00401
00402
00403
00404 #endif
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435 for (n=0; n<(NWAVE+NSCALARS); n++) {
00436 #ifdef CYLINDRICAL
00437 if (dir==1) {
00438 Wim1h[i+1][n] = (0.5*(pW[i+1][n]*r[i+1] + pW[i][n]*r[i])
00439 - (dWm[i+1][n]*r[i+1] - dWm[i][n]*r[i])/6.0)/ri[i+1];
00440 }
00441 else
00442 #endif
00443 {
00444 Wim1h[i+1][n] = 0.5*(pW[i+1][n]+pW[i][n]) - (dWm[i+1][n]-dWm[i][n])/6.0;
00445 }
00446 }
00447
00448
00449
00450
00451 for (n=0; n<(NWAVE+NSCALARS); n++) {
00452 Wlv[n] = Wim1h[i ][n];
00453 Wrv[n] = Wim1h[i+1][n];
00454 }
00455
00456
00457
00458
00459
00460 gamma_curv = 0.0;
00461 #ifdef CYLINDRICAL
00462 if (dir==1) gamma_curv = dx/(6.0*r[i]);
00463 #endif
00464
00465 for (n=0; n<(NWAVE+NSCALARS); n++) {
00466 qa = (Wrv[n]-pW[i][n])*(pW[i][n]-Wlv[n]);
00467 qb = Wrv[n]-Wlv[n];
00468 qc = 6.0*(pW[i][n] - 0.5*(Wlv[n]*(1.0-gamma_curv) + Wrv[n]*(1.0+gamma_curv)));
00469 if (qa <= 0.0) {
00470 Wlv[n] = pW[i][n];
00471 Wrv[n] = pW[i][n];
00472 } else if ((qb*qc) > (qb*qb)) {
00473 Wlv[n] = (6.0*pW[i][n] - Wrv[n]*(4.0+3.0*gamma_curv))/(2.0-3.0*gamma_curv);
00474 } else if ((qb*qc) < -(qb*qb)) {
00475 Wrv[n] = (6.0*pW[i][n] - Wlv[n]*(4.0-3.0*gamma_curv))/(2.0+3.0*gamma_curv);
00476 }
00477 }
00478
00479 for (n=0; n<(NWAVE+NSCALARS); n++) {
00480 Wlv[n] = MAX(MIN(pW[i][n],pW[i-1][n]),Wlv[n]);
00481 Wlv[n] = MIN(MAX(pW[i][n],pW[i-1][n]),Wlv[n]);
00482 Wrv[n] = MAX(MIN(pW[i][n],pW[i+1][n]),Wrv[n]);
00483 Wrv[n] = MIN(MAX(pW[i][n],pW[i+1][n]),Wrv[n]);
00484 }
00485
00486
00487
00488
00489 for (n=0; n<(NWAVE+NSCALARS); n++) {
00490 dW[n] = Wrv[n] - Wlv[n];
00491 W6[n] = 6.0*(pW[i][n] - 0.5*(Wlv[n]*(1.0-gamma_curv) + Wrv[n]*(1.0+gamma_curv)));
00492 }
00493
00494
00495
00496
00497
00498
00499 pWl = (Real *) &(Wl[i+1]);
00500 pWr = (Real *) &(Wr[i]);
00501
00502 #ifndef CTU_INTEGRATOR
00503
00504 for (n=0; n<(NWAVE+NSCALARS); n++) {
00505 pWl[n] = Wrv[n];
00506 pWr[n] = Wlv[n];
00507 }
00508
00509
00510 #elif defined(HLLE_FLUX) || defined(HLLC_FLUX) || defined(HLLD_FLUX)
00511 for (n=0; n<NWAVE; n++) {
00512 pWl[n] = Wrv[n];
00513 pWr[n] = Wlv[n];
00514 }
00515
00516 #ifdef HLL_ALL_WAVE
00517 hllallwave_flag = 1;
00518 #endif
00519
00520 for (n=0; n<NWAVE; n++) {
00521 qa = 0.0;
00522 if (hllallwave_flag || ev[n] > 0.0) {
00523 qx1 = 0.5*dtodx*ev[n];
00524 qb = qx1;
00525 qc = FOUR_3RDS*SQR(qx1);
00526 #ifdef CYLINDRICAL
00527 if (dir==1) {
00528 qxx1 = SQR(qx1)*dx/(3.0*(ri[i+1]-dx*qx1));
00529 qb -= qxx1;
00530 qc -= 2.0*qx1*qxx1;
00531 }
00532 #endif
00533 for (m=0; m<NWAVE; m++) {
00534 qa += lem[n][m]*(qb*(dW[m]-W6[m]) + qc*W6[m]);
00535 }
00536 for (m=0; m<NWAVE; m++) pWl[m] -= qa*rem[m][n];
00537 }
00538
00539 qa = 0.0;
00540 if (hllallwave_flag || ev[n] < 0.0) {
00541 qx2 = 0.5*dtodx*ev[n];
00542 qb = qx2;
00543 qc = FOUR_3RDS*SQR(qx2);
00544 #ifdef CYLINDRICAL
00545 if (dir==1) {
00546 qxx2 = SQR(qx2)*dx/(3.0*(ri[i]-dx*qx2));
00547 qb -= qxx2;
00548 qc -= 2.0*qx2*qxx2;
00549 }
00550 #endif
00551 for (m=0; m<NWAVE; m++) {
00552 qa += lem[n][m]*(qb*(dW[m]+W6[m]) + qc*W6[m]);
00553 }
00554 for (m=0; m<NWAVE; m++) pWr[m] -= qa*rem[m][n];
00555 }
00556 }
00557
00558
00559 #else
00560
00561 qx1 = 0.5*MAX(ev[NWAVE-1],0.0)*dtodx;
00562 qxx1 = 0.0;
00563 #ifdef CYLINDRICAL
00564 if (dir==1)
00565 qxx1 = SQR(qx1)*dx/(3.0*(ri[i+1]-dx*qx1));
00566 #endif
00567 for (n=0; n<(NWAVE+NSCALARS); n++) {
00568 pWl[n] = Wrv[n] - qx1 *(dW[n] - (1.0-FOUR_3RDS*qx1)*W6[n])
00569 + qxx1*(dW[n] - (1.0- 2.0*qx1)*W6[n]);
00570 }
00571
00572 qx2 = -0.5*MIN(ev[0],0.0)*dtodx;
00573 qxx2 = 0.0;
00574 #ifdef CYLINDRICAL
00575 if (dir==1)
00576 qxx2 = SQR(qx2)*dx/(3.0*(ri[i]+dx*qx2));
00577 #endif
00578 for (n=0; n<(NWAVE+NSCALARS); n++) {
00579 pWr[n] = Wlv[n] + qx2 *(dW[n] + (1.0-FOUR_3RDS*qx2)*W6[n])
00580 + qxx2*(dW[n] + (1.0- 2.0*qx2)*W6[n]);
00581 }
00582
00583
00584
00585
00586
00587
00588
00589
00590 for (n=0; n<NWAVE; n++) {
00591 if (ev[n] >= 0.0) {
00592 qa = 0.0;
00593 qx1 = 0.5*dtodx*ev[NWAVE-1];
00594 qx2 = 0.5*dtodx*ev[n];
00595 qb = qx1 - qx2;
00596 qc = FOUR_3RDS*(SQR(qx1) - SQR(qx2));
00597 #ifdef CYLINDRICAL
00598 if (dir==1) {
00599 qxx1 = SQR(qx1)*dx/(3.0*(ri[i+1]-dx*qx1));
00600 qxx2 = SQR(qx2)*dx/(3.0*(ri[i+1]-dx*qx2));
00601 qb -= qxx1 - qxx2;
00602 qc -= 2.0*(qx1*qxx1 - qx2*qxx2);
00603 }
00604 #endif
00605 for (m=0; m<NWAVE; m++) {
00606 qa += lem[n][m]*(qb*(dW[m]-W6[m]) + qc*W6[m]);
00607 }
00608 for (m=0; m<NWAVE; m++) pWl[m] += qa*rem[m][n];
00609 }
00610 }
00611
00612 for (n=0; n<NWAVE; n++) {
00613 if (ev[n] <= 0.0) {
00614 qa = 0.0;
00615 qx1 = 0.5*dtodx*ev[0];
00616 qx2 = 0.5*dtodx*ev[n];
00617 qb = qx1 - qx2;
00618 qc = FOUR_3RDS*(SQR(qx1) - SQR(qx2));
00619 #ifdef CYLINDRICAL
00620 if (dir==1) {
00621 qxx1 = SQR(qx1)*dx/(3.0*(r[i]-dx*qx1));
00622 qxx2 = SQR(qx2)*dx/(3.0*(r[i]-dx*qx2));
00623 qb -= qxx1 - qxx2;
00624 qc -= 2.0*(qx1*qxx1 - qx2*qxx2);
00625 }
00626 #endif
00627 for (m=0; m<NWAVE; m++) {
00628 qa += lem[n][m]*(qb*(dW[m]+W6[m]) + qc*W6[m]);
00629 }
00630 for (m=0; m<NWAVE; m++) pWr[m] += qa*rem[m][n];
00631 }
00632 }
00633
00634
00635 for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00636 if (W[i].Vx > 0.) {
00637 qb = 0.5*dtodx*(ev[NWAVE-1]-W[i].Vx);
00638 qc = 0.5*dtodx*dtodx*TWO_3RDS*(SQR(ev[NWAVE-1]) - SQR(W[i].Vx));
00639 pWl[n] += qb*(dW[m]-W6[m]) + qc*W6[m];
00640 } else if (W[i].Vx < 0.) {
00641 qb = 0.5*dtodx*(ev[0]-W[i].Vx);
00642 qc = 0.5*dtodx*dtodx*TWO_3RDS*(ev[0]*ev[0] - W[i].Vx*W[i].Vx);
00643 pWr[n] += qb*(dW[m]+W6[m]) + qc*W6[m];
00644 }
00645 }
00646
00647 #endif
00648
00649
00650
00651
00652 for (m=0; m<NWAVE; m++) {
00653 ev[m] = ev_ip1[m];
00654 for (n=0; n<NWAVE; n++) {
00655 rem[m][n] = rem_ip1[m][n];
00656 lem[m][n] = lem_ip1[m][n];
00657 }
00658 }
00659
00660 }
00661
00662 return;
00663 }
00664
00665
00666
00667
00668
00669
00670 void lr_states_init(MeshS *pM)
00671 {
00672 int nmax,size1=0,size2=0,size3=0,nl,nd;
00673
00674
00675 for (nl=0; nl<(pM->NLevels); nl++){
00676 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00677 if (pM->Domain[nl][nd].Grid != NULL) {
00678 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00679 size1 = pM->Domain[nl][nd].Grid->Nx[0];
00680 }
00681 if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
00682 size2 = pM->Domain[nl][nd].Grid->Nx[1];
00683 }
00684 if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
00685 size3 = pM->Domain[nl][nd].Grid->Nx[2];
00686 }
00687 }
00688 }
00689 }
00690
00691 size1 = size1 + 2*nghost;
00692 size2 = size2 + 2*nghost;
00693 size3 = size3 + 2*nghost;
00694 nmax = MAX((MAX(size1,size2)),size3);
00695
00696 if ((pW = (Real**)malloc(nmax*sizeof(Real*))) == NULL) goto on_error;
00697
00698 if ((dWm = (Real**)calloc_2d_array(nmax, NWAVE, sizeof(Real))) == NULL)
00699 goto on_error;
00700
00701 if ((Wim1h = (Real**)calloc_2d_array(nmax, NWAVE, sizeof(Real))) == NULL)
00702 goto on_error;
00703
00704 return;
00705 on_error:
00706 lr_states_destruct();
00707 ath_error("[lr_states_init]: malloc returned a NULL pointer\n");
00708 }
00709
00710
00711
00712
00713
00714 void lr_states_destruct(void)
00715 {
00716 if (pW != NULL) free(pW);
00717 if (dWm != NULL) free_2d_array(dWm);
00718 if (Wim1h != NULL) free_2d_array(Wim1h);
00719 return;
00720 }
00721
00722 #endif