| 162 |
|
static float * |
| 163 |
|
getambdiffs(AMBHEMI *hp) |
| 164 |
|
{ |
| 165 |
< |
const double normf = 1./bright(hp->acoef); |
| 165 |
> |
const double normf = 1./(pbright(hp->acoef) + FTINY); |
| 166 |
|
float *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float)); |
| 167 |
< |
float *ep; |
| 167 |
> |
float *ep, *earr2; |
| 168 |
|
AMBSAMP *ap; |
| 169 |
|
double b, b1, d2; |
| 170 |
|
int i, j; |
| 198 |
|
ep[-hp->ns-1] += d2; |
| 199 |
|
} |
| 200 |
|
/* correct for number of neighbors */ |
| 201 |
< |
earr[0] *= 8./3.; |
| 202 |
< |
earr[hp->ns-1] *= 8./3.; |
| 203 |
< |
earr[(hp->ns-1)*hp->ns] *= 8./3.; |
| 204 |
< |
earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 8./3.; |
| 201 |
> |
earr[0] *= 6./3.; |
| 202 |
> |
earr[hp->ns-1] *= 6./3.; |
| 203 |
> |
earr[(hp->ns-1)*hp->ns] *= 6./3.; |
| 204 |
> |
earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 6./3.; |
| 205 |
|
for (i = 1; i < hp->ns-1; i++) { |
| 206 |
< |
earr[i*hp->ns] *= 8./5.; |
| 207 |
< |
earr[i*hp->ns + hp->ns-1] *= 8./5.; |
| 206 |
> |
earr[i*hp->ns] *= 6./5.; |
| 207 |
> |
earr[i*hp->ns + hp->ns-1] *= 6./5.; |
| 208 |
|
} |
| 209 |
|
for (j = 1; j < hp->ns-1; j++) { |
| 210 |
< |
earr[j] *= 8./5.; |
| 211 |
< |
earr[(hp->ns-1)*hp->ns + j] *= 8./5.; |
| 210 |
> |
earr[j] *= 6./5.; |
| 211 |
> |
earr[(hp->ns-1)*hp->ns + j] *= 6./5.; |
| 212 |
|
} |
| 213 |
+ |
/* preen map to avoid cliffs */ |
| 214 |
+ |
earr2 = (float *)malloc(hp->ns*hp->ns*sizeof(float)); |
| 215 |
+ |
if (earr2 == NULL) |
| 216 |
+ |
return(earr); |
| 217 |
+ |
memcpy(earr2, earr, hp->ns*hp->ns*sizeof(float)); |
| 218 |
+ |
for (i = 0; i < hp->ns-1; i++) { |
| 219 |
+ |
float *ep2 = earr2 + i*hp->ns; |
| 220 |
+ |
ep = earr + i*hp->ns; |
| 221 |
+ |
for (j = 0; j < hp->ns-1; j++, ep2++, ep++) { |
| 222 |
+ |
if (ep2[1] < .5*ep2[0]) { |
| 223 |
+ |
ep[0] -= .125*ep2[0]; |
| 224 |
+ |
ep[1] += .125*ep2[0]; |
| 225 |
+ |
} else if (ep2[1] > 2.*ep2[0]) { |
| 226 |
+ |
ep[1] -= .125*ep2[1]; |
| 227 |
+ |
ep[0] += .125*ep2[1]; |
| 228 |
+ |
} |
| 229 |
+ |
if (ep2[hp->ns] < .5*ep2[0]) { |
| 230 |
+ |
ep[0] -= .125*ep2[0]; |
| 231 |
+ |
ep[hp->ns] += .125*ep2[0]; |
| 232 |
+ |
} else if (ep2[hp->ns] > 2.*ep2[0]) { |
| 233 |
+ |
ep[hp->ns] -= .125*ep2[hp->ns]; |
| 234 |
+ |
ep[0] += .125*ep2[hp->ns]; |
| 235 |
+ |
} |
| 236 |
+ |
} |
| 237 |
+ |
} |
| 238 |
+ |
free(earr2); |
| 239 |
|
return(earr); |
| 240 |
|
} |
| 241 |
|
|
