Array of functions and Segmentation fault - invalid memory reference

I am trying to set my function f as an array, but I get the following error:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:

#0  0x6f8b36e3
#1  0x6f8a2722
#2  0x402752
#3  0x747bd411

I have to solve the Kepler's equation:f=psi-e*sin(psi)-M for each value of M.So, if I have an array M of dimension 8, my program will calculate 8 zeros. The thing is that, if I write f=psi-e*sin(psi)-M(1) I will calculate the first zero,and if I write f=psi-e*sin(psi)-M(8) I will calculate the last zero.But ,my problem is that if I want to calculate all zeros at once, I would have to write f=psi-e*sin(psi)-M(1:8) and my program should type all zeros, but,that doesnt happen and i get the error i mentioned earlier.Here is the code:

SUBROUTINE (I USED EXTERNALLY):This subroutine is the bisection method (to get zeros):

   subroutine bisecc(f,xl,xr,kmax,tol,k,xm)
  implicit real*8 (a-h,o-z)
  real*8 f
  fl=f(xl)
  fr=f(xr)

  if(fl*fr .gt. 0.0D0) goto 100  
  do k=1,kmax
    xm=(xr+xl)/2.0D0              
    fm=f(xm)

    dif=abs((xr-xl)/xm)        
    if(dif .lt. tol) goto 200    
    write(*,*) k,xm!,dif  
    if (fm*fr .le. 0.0D0) then
      xl=xm
      fl=fm
    else
      xr=xm
      fr=fm
    end if
  end do

  return
200   write(*,*) 'WISHED PRECISION  REACHED'
      return
100   write(*,*) 'BAD CHOICE OF DATA'
      return
      end

MAIN PROGRAM:

       include 'bisecc.f' 
   implicit real*8 (a-h,o-z)            
   external f     
   real*8 f    
   ! I  WRITE  THE INTERVAL OF MY 8 ZEROS(left and right point)
   b=0.1D0


   xl1=-0.5D0                
   xr1=0.D0                 
   xl2=xr1+b
   xr2=1.D0
   xl3=xr2+b
   xr3=2.D0
   xl4=xr3+b
   xr4=3.D0
   xl5=xr4+b
   xr5=4.D0
   xl6=xr5+b
   xr6=5.D0
   xl7=xr6+b
   xr7=6.D0
   xl8=xr7+b
   xr8=7.D0        
   kmax=100                            
   tol=0.0001D0                               
   call bisecc(f,xl1,xr1,kmax,tol,k,xm1)         
   call bisecc(f,xl2,xr2,kmax,tol,k,xm2)          
   call bisecc(f,xl3,xr3,kmax,tol,k,xm3)         
   call bisecc(f,xl4,xr4,kmax,tol,k,xm4)         
   call bisecc(f,xl5,xr5,kmax,tol,k,xm5)         
   call bisecc(f,xl6,xr6,kmax,tol,k,xm6)         
   call bisecc(f,xl7,xr7,kmax,tol,k,xm7)          
   call bisecc(f,xl8,xr8,kmax,tol,k,xm8)                     
   write(*,*) 'Program ended'
   stop
   end program

   real*8 function f(psi)
   implicit real*8 (a-h,o-z)
   real*8 M(8)
   dimension f(8)
      e=0.2056D0                          
      pi=acos(-1.0D0)                        
      M=(/pi/4.D0,pi/2.D0,3.D0/4.D0*pi,pi,5.D0/4.D0*pi,3.D0*
 &    pi/2.D0,7.D0/4.D0*pi,2.D0*pi/)
      c=sqrt((1.0D0-e)/(1.0D0+e))                
   f=psi-e*sin(psi)-M(1:8)             !KEPLER EQUATION

   return
   end function

EXAMPLE: Here I wanted to calculate the value of psi for the first value of M, M(1)=pi/4.

In https://i.sstatic.net/qxEv2.jpg you can see that psi=0.95303344726562489. So I have just calculated the first zero. But, you can also see this message 7 times datos mal elegidos. It means that the program can only show me that zero (for M(1)), and the other 7 zeros are not calculated, because I wrote f=psi-e*sin(psi)-M(1). What should I write so I can get the result of all zeros instead of 1, like in this example?

Solution

Because the function f() is used in the bisection routine bisecc(), I think it would be much simpler to pass each input to bisecc() via a DO loop, rather than making f() a function returning an array (because the latter requires to modify bisecc() also). We can pass the value of M to f() in various ways (which is almost FAQ and I believe there are a lot of Q/A pages). One simple way is to contain f() in the main program and use host association for M. So a simplified code may look like

program main
    implicit none
    integer  kmax, kiter, i
    real*8   xl( 8 ), xr( 8 ), xans( 8 ), tol, M( 8 ), b, pi

    pi   = acos(-1.0D0)                 
    kmax = 100
    tol  = 1.0d-8

    M = [ pi/4.D0,      pi/2.D0,      3.D0/4.D0*pi, pi, &
          5.D0/4.D0*pi, 3.D0*pi/2.D0, 7.D0/4.D0*pi, 2.D0*pi ]
    ! or M = [( i, i=1,8 )] * pi/4.0D0

    ! Use a fixed interval for simplicity.
    xl   =  0.0d0
    xr   = 10.0d0
    xans =  0.0d0

    do i = 1, 8
        call bisecc( f, xl( i ), xr( i ), kmax, tol, kiter, xans( i ) )
        ! print *, "check: f(xans(i)) = ", f( xans( i ) )
    enddo

contains

function f( psi ) result( res )
    implicit none
    real*8  psi, e, res
    e = 0.2056D0
    res = psi - e * sin( psi ) - M( i )   !<-- this "M(i)" refers to that defined above
end function 

end program

with an external bisecc routine (a little modified so as not to use GOTO)

subroutine bisecc( f, xl, xr, kmax, tol, k, xm )
    implicit none
    real*8  f, xl, xr, tol, xm
    external f
    integer kmax, k
    real*8  fl, fr, fm, dif

    fl = f( xl )
    fr = f( xr )
    if( fl * fr > 0.0D0 ) then
        write(*,*) "bad input data (xl,xr)"
        return
    endif

    do k = 1, kmax
        xm = (xr + xl) / 2.0D0              
        fm = f( xm )

        dif = abs( (xr-xl) / xm )
        if ( dif < tol ) then
            write(*,*) "bisection converged: k=", k, "xm=", xm
            return
        endif

        if ( fm * fr <= 0.0D0 ) then
            xl = xm
            fl = fm
        else
            xr = xm
            fr = fm
        end if
    end do  !! iteration

    write(*,*) "bisection did not converge: k=", k, "xm=", xm
end

which gives

 bisection converged: k=          31 xm=  0.95299366395920515     
 bisection converged: k=          31 xm=   1.7722388869151473     
 bisection converged: k=          30 xm=   2.4821592587977648     
 bisection converged: k=          30 xm=   3.1415926571935415     
 bisection converged: k=          29 xm=   3.8010260276496410     
 bisection converged: k=          29 xm=   4.5109464414417744     
 bisection converged: k=          29 xm=   5.3301916457712650     
 bisection converged: k=          29 xm=   6.2831853143870831

The answer seems to agree with the plot of the Kepler equation with e = 0.2056 (so bisecc() is probably OK).

The above code still has a lot of points for improvement. In particular, it is usually more convenient to include a function like f() into a module (or even include all routines into a module). We can also pass M by making it a module variable and use it from f() (rather than using common statements) or via host association, so please try it if interested.