      implicit none
      real*8 u2,u3,u4,u5
      real*8 y,miny,maxy,dy,jqy,jia,j1,j2,j3
      character*50 jtfile,jiafile,j1file,j2file,j3file

      write (*,'(T5,A,$)') 'Enter min, max and interval of Y ::'
      read (*,*) miny,maxy,dy
      write (*,'(T5,A,$)') 'Enter a2 (1/A^2), b3 (1/A^3), a4 (1/A^4) and b5 (1/A^5):'
      read (*,*) u2,u3,u4,u5


      jtfile = "jaa.dat"
      jiafile = "jia.dat"
      j1file = "j1.dat"
      j2file = "j2.dat"
      j3file = "j3.dat"

      open (unit=10,file=jtfile,status='unknown')
      open (unit=11,file=jiafile,status='unknown')
      open (unit=12,file=j1file,status='unknown')
      open (unit=13,file=j2file,status='unknown')

      INQUIRE (unit=10,name=jtfile)
      INQUIRE (unit=11,name=jiafile)
      INQUIRE (unit=12,name=j1file)
      INQUIRE (unit=13,name=j2file)

      if (u5 .gt. 0.0) then 
         open (unit=14,file=j3file,status='unknown')
         INQUIRE (unit=14,name=j3file)
      end if

c
c Calculate AA components and save to file
c

      y = miny
      do while (y .LT. maxy)
         CALL jaa_cal(y,u2,u3,u4,u5,jia,j1,j2,j3)
         jqy = jia+j1+j2+j3
         write (10,'(3G11.4)') y,jqy,0.05*jqy
         write (11,'(3G11.4)') y,jia,0.05*jia
         write (12,'(3G11.4)') y,j1,0.05*j1
         write (13,'(3G11.4)') y,j2,0.05*j2
         if (u5 .gt. 0.0) write (14,'(3G11.4)') y,j3,0.05*j3
         y = y + dy
      end do

	write (*,*) ' '
	write (*,*) 'Sum stored in '//jtfile
	write (*,*) 'Jia stored in '//jiafile
	write (*,*) 'J1 stored in '//j1file
	write (*,*) 'J2 stored in '//j2file
	if (u5 .gt. 0.0) write (*,*) 'J3 stored in '//j3file
	STOP
	END



      SUBROUTINE jaa_cal(y,u2,u3,u4,u5,jia,j1,j2,j3)
c
c Additive Approach Model function (in the incoherent limit).
c
      implicit none
      real*8 y,wd,u2,u3,u4,u5
      real*8 jia,j1,j2,j3,pi

      pi = 4.0*ATAN(1.0)

      wd = Y/SQRT(u2)
      Jia= 1/SQRT(2*pi*u2)*DEXP(-Y**2/(2*u2))
      J1 = -u3*wd/(2*U2**1.5)*(1 - wd**2/3)*Jia
      J2 = u4/(8*u2**2)*(1 - 2*wd**2 + wd**4/3)*Jia
      J3 = u5*wd/(8*u2**2.5)*(1 - 2/3*wd**2 + wd**4/15)*Jia

      return
      end
