Opened 11 years ago

Last modified 11 years ago

#18 new defect - etc

code keeps track of exp(-tau) incorrectly?

Reported by: peter Owned by: nobody
Priority: critical Milestone:
Component: etc Version:
Keywords: Cc:


When one adds the following print statement to e2() in service.cpp

if( fabs(exp(-t)-exp_mt)/SDIV(exp_mt) > 0.001 )

printf( "COMP e2: %.16e %.16e %.16e\

", t, exp(-t), exp_mt );

one gets the following output (excerpt):

COMP e2: 1.9192649841308594e+01 4.6210224754419657e-09 8.9803414704547890e-11
COMP e2: 1.5586171150207520e+01 1.7022050533117493e-07 6.9368506494527082e-09
COMP e2: 1.2656538009643555e+01 3.1866574569759630e-06 2.3698864026755473e-07
COMP e2: 1.0276868820190430e+01 3.4420135855679434e-05 4.1725779738044366e-06
COMP e2: 8.3440465927124023e+00 2.3780808081795259e-04 4.2871797631960362e-05
COMP e2: 6.7742619514465332e+00 1.1428136407263399e-03 2.8438138542696834e-04
COMP e2: 5.4994125366210938e+00 4.0891729723596662e-03 1.3220318360254169e-03
COMP e2: 4.4641528129577637e+00 1.1514446518938374e-02 4.6043256297707558e-03
COMP e2: 3.6235115528106689e+00 2.6688792645240801e-02 1.2682755477726460e-02
COMP e2: 2.9409501552581787e+00 5.2815521944114145e-02 2.8874300420284271e-02
COMP e2: 2.3867785930633545e+00 9.1925336329503518e-02 5.6312073022127151e-02
COMP e2: 1.9378761053085327e+00 1.4400948618249837e-01 9.6735745668411255e-02
COMP e2: 1.5762796401977539e+00 2.0674282690849677e-01 1.4957897365093231e-01
COMP e2: 1.2820999622344971e+00 2.7745404524375594e-01 2.1323759853839874e-01
COMP e2: 1.0419369935989380e+00 3.5277070514551945e-01 2.8482690453529358e-01

there are clearly significant differences between exp(-t) and exp_mt, which should in fact be nearly identical. This could imply that the code keeps track of exp(-tau) incorrectly, or that the e2 is called with the wrong second argument. Either way it would be an important bug.

Change History (3)

comment:1 Changed 11 years ago by peter

Type: defect - FPEdefect - etc

comment:2 Changed 11 years ago by peter

PS - the above output was obtained by running

comment:3 Changed 11 years ago by peter

The origin of this problems was in calls of the type:

E2 = e2(opac.TauAbsFace[i],opac.ExpmTau[i]);

where opac.ExpmTau[i] != exp(-opac.TauAbsFace[i]). One problem is that opac.ExpmTau[i] contains an extra factor geometry.DirectionalCosin which is absent in opac.TauAbsFace[i]. Hence the second parameter to e2 has been removed in r793. But there is a second problem, the origin of which seems to be in the routine IterStart() in iter_startend.cpp. This contains the following code in r799:

opac.TauAbsFace[i] = opac.taumin;
opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin);
opac.ExpmTau[i] = (float)opac.ExpZone[i];

Here opac.opacity_abs[i]*radius.drad_x_fillfac/2. can be vastly different from opac.taumin even when geometry.DirectionalCosin == 1. (which is usually the case).

Note: See TracTickets for help on using tickets.