What's going on here? Vector does not want to receive a number

Help
1

So here’s the thing, for some reason if i calculate the value for an array component, it just don’t recieve the value and return it as NaN, but if i calculate the same value and define a variable with this number, i can print this number and it is right, but again, if i put this number in the vector, it recieve NaN, i have no idea what’s going on. I need to calculate g(2) ang g(3), if i do the same calculation but don’t put it in the array, the number is correct, but it just doesn’t want to enter the array

           grl(1) = g(1)*cosx + sqrt(g(2)**2 + g(3)**2)*sin(e)*sinx
           
           grl(2) = gr*sinx*sin(e) !old method that works but need to change
           !g2 = g(2)*cosx + ((gr*g(3)*cos(e) - g(1)*g(2)*sin(e))*sinx) / sqrt(g(2)**2 + g(3)**2) !the correct way to calculate, but in a variable because the array doesnt want to receive it
            
           grl(3) = gr*sinx*cos(e) !old method that works but need to change
           !g3 = g(3)*cosx - ((gr*g(2)*cos(e) - g(1)*g(3)*sin(e))/(sqrt(g(2)**2 + g(3)**2)))*sinx !the same calculation but in a variable
            
           !print*,g(2)*cosx + ((gr*g(3)*cos(e) - g(1)*g(2)*sin(e))*sinx) / sqrt(g(2)**2 + g(3)**2) !if i print it the number is correct
           !print*, g2 !print the correct number
           grl(2) = g2 !if i put the correct number in the array IT RECIEVES NAN FOR SOME REASON
           print*, grl(2)` "NAN NAN NAN ..."

The entire code:
hs.f90 (7.0 KB)
mmc1-he3.csv (1.3 MB)

2

Can you attach the entire source code or show the variable definitions?

g2 seems to be commented so… isn’t the case that you are forgetting to initialize it or maybe the values in your calculation are wrong? Check the variable “sinx” or the values in the vector “g”

I can’t help much if you don’t show more context code.

3

I notice that you have this in your expressions. If you are generating NaN values or overflows or underflows, you might consider replacing this with hypot(g(2),g(3)). That latter is an intrinsic function (f2008 and later) that is supposed to avoid problems with overflows (for large values) and underflows (for small values). If you are using a compiler that does not have that intrinsic, then a google search will show a few replacements written in fortran (e.g. there is one in the LAPACK library). The correct way to compute that value is a little complicated, involving prescaling, some tests for special cases, the computation, and then post scaling.

To better answer your full question, you should show more of the code, including the declarations.

1 Like
4

sure, here it is
hs.f90 (7.0 KB)
mmc1-he3.csv (1.3 MB)

5

So… as I suggested to you before keep compiling using debug flags I mentioned while you are developing, like:

gfortran -g -fcheck=bounds -Wall -Wextra -fmax-errors=1 -fbacktrace

As you can see the problem is in another castle:

                                           NaN                                           NaN
                                           NaN                                           NaN
                                           NaN                                           NaN
                                           NaN                                           NaN
                                           NaN                                           NaN
At line 151 of file /home/ian/downloads/hs (1).f90
Fortran runtime error: Index '0' of dimension 1 of array 'sigma' below lower bound of 1

Error termination. Backtrace:
#0  0x55bff1284344 in rgd2
	at /home/ian/downloads/hs (1).f90:151
#1  0x55bff1285b66 in main
	at /home/ian/downloads/hs (1).f90:221

The big problem is the calculation of matrixj reaching zero: as you have a log(1.0 + x) in your expression, x can’t be 0 as \log(1) = 0 and as you will see for yourself it is.

What happens when both particles have the same velocity, do your calculation holds?

1 Like
6

but the code works fine when i use the old method to calculate the post collisions velocities, using the matrixj as it is, so i don’t understand

7

Your calculation uses data generated using a pseudo-random number generator. Therefore, errors may or may not occur in an unpredictable manner in repeated runs with no change to the code at all.

Your program contains the following lines:

            !CALCULATING THE INDEX J OF THE MATRIX
            matrixj =  FLOOR(log(1d0 + (gr*v0)/400d0)/(log(1.005)) + 1d0/2d0)
           ! write(*,*)'j = ', matrixj
            if(matrixj > 800)matrixj = 800

You should probably add this:

if (matrixj < 1) matrixj = 1

to avoid using a zero subscript (whether that happens depends on the RNG that you use).

In general, you should turn subscript checking on, and use other error checking facilities available with your compiler until you are confident that your code is error free.

1 Like
8

Your new method is wrong, it is not “that the vector doesn’t want to receive a number”, your calculation of g2 relies on an expression in the denominator, sqrt(g(2)**2 + g(3)**2), if that thing gets to zero you’ll get a NaN in g2.

Funnily enough, you are getting the components g(2) and g(3) that are exactly zero:

g  =   -2.00000000000000000000000000000000000         0.00000000000000000000000000000000000         0.00000000000000000000000000000000000      
g2 =                                            NaN

Is that part supposed to be the norm of the vector g? if so you can use the function norm2

An ancient debug method is printing all the variables involved

print*, 'sqrt=', sqrt(g(2)**2 + g(3)**2)
print*, 'x  = ', (gr*v0)/400d0
print*, 'gr = ', gr
print*, 'g  = ', g
print*, 'g2 = ', g2
print*, 'v0 = ', v0

Hope it helps now.

1 Like
9

The Taylor expansion for log(1+x) is x-x**2/2+... Evaluating it as log(1.0+x) will give no correct significant digits when x is smaller than epsilon, whereas evaluating it with the low order Taylor expansion will give all correct significant digits. I can’t tell from the discussion here if that is a problem, but it might be worth looking at. Some numerical libraries have separate functions for log(x) and log(1+x), but standard fortran does not provide the latter function, so you will need to write your own if necessary.

1 Like
11

Wow, how didn’t i see that. Thank you!! i added a line if(g(2)==0.0 .and. g(3)==0.0)cycle and now it seems to work, but i’ll have to ask my advisor tomorrow about what happens to the sqrt(g(2)**2 + g(3)**2) term when g(2) and g(3) are both 0.

2 Likes