ERROR: LoadError: DomainError: in quaternion(::Array{Float64,2})

I get the error message

ERROR: LoadError: DomainError: in quaternion(::Array{Float64,2})

when I run the the following script (test5.jl):

function quaternion(m)

       q = [0.0,0.0,0.0,0.0]                   
         
       q[1] = 0.5*sqrt(1+ m[1,1] + m[2,2]+ m[3,3])                
        if q[1] == 0.0                                            
         q[1] = eps()                                                 
        end 
       q[2] = (1/(4*q[1]))*(m[3,2] - m[2,3])
       q[3] = (1/(4*q[1]))*(m[1,3] - m[3,1])
       q[4] = (1/(4*q[1]))*(m[2,1] - m[1,2])

       return q

end

lat = 8.7758905947561e-8
lon = 0.0

C_ne(lat,lon) 	= [-8.77589e-8 -0.0 -1.0; -0.0 1.0 -0.0; 1.0 0.0 -8.77589e-8]
C__eb		= [0.5 0.0 -1.0; 0.0 1.0 0.0; 1.0 0.0 2.22045e-16]

C__nb = C__eb * C_ne(lat,lon)
# C__nb= [-1.0 0.0 -0.5; 0.0 1.0 -0.0; -8.77589e-8 0.0 -1.0]

s = quaternion(C__nb)

println("s= ",s)

But if I switch to (i.e. load C__nb manually)

<code>

# C__nb = C__eb * C_ne(lat,lon)
C__nb= [-1.0 0.0 -0.5; 0.0 1.0 -0.0; -8.77589e-8 0.0 -1.0]

<code>

The code executes fine:

$ julia test5.jl 
s= [2.22045e-16,0.0,-5.6295e14,0.0]

Why do I get an execution error in my function quaternion(m) when the matrix C__nb is computed?

Because you are taking the sqrt of a negative real number. Convert it to a complex number first if it can be negative.

OK now I see,

When I compute C__nb the value of (1+ m[1,1] + m[2,2]+ m[3,3]) becomes: -4.38794500734474e-8

But when I copy the result of the expression C__nb = C__eb * C_ne(lat,lon) from a println(C__nb) and use it “manually”, the evaluation of (1+ m[1,1] + m[2,2]+ m[3,3]) becomes: 0.0

Suttle…

Thanks!