-----------------------
     INTRODUCTION
-----------------------
                                   FLUX.F
 		             Laura E. Marcucci
                     University of Pisa and INFN - sezione di Pisa

   Code for the calculation of the proton-proton weak capture reaction
   rate [cm^3 mol^(-1) s^(-1)] presented in Tognelli, Degl'Innocenti,
   Marcucci, Prada Moroni, 2015, Physics Letters B 742 (2015), pp. 189-194.

   Two subroutines (flux and b5) are needed to calculate the flux for the
   proton-proton weak capture reaction, given the star temperature in
   10^9 K. It is based on the following articles:

    - L.E. Marcucci et al., Phys. Rev. Lett. 110, 192503 (2013) 

    - NACRE II, Nucl. Phys. A 918, 61 (2013) [ Eq.(3) ].


   The code is structured in such a way that there are several options:

   (1) the initial p-p wave function can include only the S-wave contribution
       or also all the P-wave contributions. The S+P waves is what has been
       suggested in L.E. Marcucci et al., the S-wave contribution is what
       has been used so far in the literature.

   (2) The user can calculate the flux using the calculated values for
       the p-p astrophysical S-factor on an energy grid of 101 points
       of step = 1 keV, or alternatively use the values of S(0), S^n(0),
       [S^n(0) is the n-th derivative of S(E) calculated at zero
       energy] up to n=3, and tabulate S(E) with the preferred steps.
       In this second case, in parameter, the user needs to adjust
       he and ne_int (see below).

-----------------------
HOW TO COMPILE THE CODE
-----------------------
   The code should be compiled as:

     (gfortran compiler)
     gfortran -O4 -o flux flux.f
     
     or

     (intel compiler)
     ifort -O4 -o flux flux.f

-----------------------
        USAGE
-----------------------
   The user needs to know only the following inputs indices:

     iprint = 1/0 for printing/not printing intermediate results

     t9 = star temperature in 10^9 K

     ic = 1 for S+P waves
        = 2 for only S-waves

     inter = 1 for using the calculated S-factor
           = 2 for using the fitted values 
               for S(0), S'(0), S''(0) and S'''(0) and then 
               tabulating the S(E) as preferred.
               In this case, in parameter the user should fix:
                         he=1.e-3 ne_int=101 : step (given in MeV) of 1 keV
                         he=1.e-4 ne_int=1001: step (given in MeV0 of 0.1 keV 
                           etc....
     
     In output, the user finds:

     flx = flux in cm^3 mol^{-1} s^{-1}
     dflx= theoretical uncertainty on the flux

-----------------------
       EXAMPLE
-----------------------

      Program example
      implicit real*8(a-h,o-z)

      t9=1.55e-2 ! sun temperature as an example (in units of 10^9 K)

      flx =0.d0
      dflx=0.d0

 1    write(*,*)'give iprint, ic, inter [see comments in flux.f]'
      write(*,*)'ic=0 to stop'

      read(*,*)iprint,ic,inter

      flx =0.d0
      dflx=0.d0

      if(ic.ne.0)then

         call flux(iprint,t9,ic,inter,flx,dflx)
         write(*,1000)flx,dflx
         go to 1

      endif

 1000 format('flx = ',e15.8,' +/- ',e15.8,/)

      end

-----------------------
         TEST
-----------------------

Running the example program one should get the following results:

give iprint, ic, inter [see comments in flux.f]
 ic=0 to stop
0 1 1
flx =  0.46154623E-19 +/-  0.61429411E-22

 give iprint, ic, inter [see comments in flux.f]
 ic=0 to stop
0 1 2
flx =  0.46263927E-19 +/-  0.72186794E-22

 give iprint, ic, inter [see comments in flux.f]
 ic=0 to stop
0 2 1   
flx =  0.45689162E-19 +/-  0.56787487E-22

 give iprint, ic, inter [see comments in flux.f]
 ic=0 to stop
0 2 2
flx =  0.45842303E-19 +/-  0.60834237E-22

 give iprint, ic, inter [see comments in flux.f]
 ic=0 to stop
0 0 0

--------------------------------------------------------------------------
     For any assistance, send an email to laura.marcucci@df.unipi.it
--------------------------------------------------------------------------