I don't understand the output of a function in a Fortran95 code

Question

I just got started with fotran95; I was given a code and I am studying it; I came across a subroutine that calls a function but I don't understand what the output is:

here is the subroutine:

SUBROUTINE collisione(vga, ga, vgb, gb)

IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: ga, gb    
DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: vga, vgb  

DOUBLE PRECISION, DIMENSION(3)   :: r, ra, rb, r_r

ra = pos_ini(ga)
rb = pos_ini(gb)


END SUBROUTINE

here is the function:

FUNCTION pos_ini(g)

IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3) :: pos_ini
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: g

DOUBLE PRECISION, DIMENSION(3) :: es, eg, es_p
DOUBLE PRECISION :: rnd, b, kos, ang 

CALL RANDOM_NUMBER(rnd)
b = 1.D0 - 2.D0*rnd 
kos = SQRT(1.D0 - b*b)
CALL RANDOM_NUMBER(rnd)
ang = pi2 * rnd

es(1) = kos * COS(ang)
es(2) = kos * SIN(ang)
es(3) = b

eg = g / SQRT(DOT_PRODUCT(g,g))

es_p = es - DOT_PRODUCT(es,eg) * eg

pos_ini = d * es_p/SQRT(DOT_PRODUCT(es_p,es_p)) ! IS THIS THE OUTPUT THAT GIVES THE VALUE OF ra and rb???????

END FUNCTION

(read the comment). So here comes my question: in the subroutine I see that the variables ra and rb are defined after using the function pos_ini where the input is the vector ga or gb. In the function pos_ini however I don't understand what the output is; is it pos_ini = d * es_p/SQRT(DOT_PRODUCT(es_p,es_p)) ?? if so why? and how come there is no intent(OUT) in function pos_ini?

Answer 1

pos_ini is a double precision function returning an array of three elements. I assume that it is declared inside a module (because there is no declaration for pi2 and for d), but that does not really matter here.

In Fortran, the return value is set via assignment to the function name (or to the name in a RESULT clause, which does not apply here). This is done with the line

pos_ini = d * es_p/SQRT(DOT_PRODUCT(es_p,es_p))

which puzzled you.

In this case, this is an array assignment. In Fortran, you can have assignments to arrays from scalars or other arrays, array arithmetic (which is done element-wise) and mixed arithmetic between scalars and arrays.

Because d lacks a definition, I cannot say what it is. It would be valid either as an array with an extent of 3 or as a scalar. es_p is an array with the same extent as pos_ini, so that is OK. DOT_PRODUCT should be self-explanatory; it is a sclar, as is SQRT.

I would recommend you to read this article on Fortran 95 features.

Answer 2

So this line

FUNCTION pos_ini(g)

begins the definition of a function called pos_ini. And this line

DOUBLE PRECISION, DIMENSION(3) :: pos_ini

tells us (and the compiler) that pos_ini will return a double precision array with 3 elements. Unless coded otherwise (with a result clause) a Fortran function declares (automatically or, as here, by explicitly writing code) a result variable of the same name as the function. And the line

pos_ini = d * es_p/SQRT(DOT_PRODUCT(es_p,es_p))

which defines a value for pos_ini defines what will be returned from an invocation of the function. There is no need, and it would be an error that the compiler would catch, to declare the return variable to have intent(out).

For functions which return a scalar value of an intrinsic type the declaration of the function's type (or of the type of the return variable if you want to think of things in that way) you can write a function opening statement along the lines of

DOUBLE PRECISION FUNCTION pos_ini(g)

but, since your function returns an array you have to code as you have, declaring the return variable explicitly.

Answer 3

In Fortran the returned function value is connected to the function name or a variable in the optional result clause. Therefore the function value will be whatever happened to be the value of pos_ini at the end of the function.

In your case it is the espression d * es_p/SQRT(DOT_PRODUCT(es_p,es_p)).

This is covered by any Fortran tutorial.