Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
Possible crippling bug in NR_(SIMP,ROMB)
I get erroneous results when using the NR_SIMP and NR_ROMB routines recursively that I belive traces to the routines themselves. Here is one of the simplest routines I could think of to illustrate the problem (the idea for this recursive calling comes from NR, section 4.6): ========= FUNCTION myfunc, x return, 2.0*x END FUNCTION innerint, y return, nr_qromb('myfunc', 0.0, 5.0) END PRO simpleint print, nr_qromb('innerint', 0.0, 2.0) END ========= Now, executing a statement like 'print, nr_qromb('myfunc',0.0,5.0)' correctly returns 25.000. However when adding the extra integral from 0-2 over y as shown above in 'simpleint', IDL returns 37.500 when the correct answer should be 50. I have run this exact calculation through the C versions of the NR routines and they correctly report 50.000. I have checked this on SGI's, Alpha's, and DecStations with the same results. It seems the recursive calling is mangling the results. (Strangely, using the y range of 0-3 gives 50!) If I replace the return statment in 'innerint' with 'return, 25.0' I (of course) get 50. As listed above, 'innerint' should be returning 25.000 every call. Of course the integral above doesn't need to be calculated this way, but I have some that I need to do using the recursive trick. Has anyone encountered this before and found a workaround? Am I missing something obvious about recursive function calling in IDL? Thanks. mh -- Matt Haffner mhaffner@fosters.astro.wisc.edu Dept. of Astronomy University of Wisconsin, Madison