      program iexp
      implicit double precision (A-H,O-Z)
      parameter (nb=11)
      dimension x(nb),y(nb),z(nb),aa(3),bb(3),sigma(nb)
      real*4 xs(nb),ys(nb),xs1(100),ys1(100),ys2(100),xstep,ystep
      real*4 as,bs,ps
      character*1 dummy
      write(*,*) 'zadaj konstanty a, b'
      read(*,*) a,b
c      a=4.d0
c      b=0.3d0
      write(*,*) 'zadaj oblast integrovania alfa, beta'
      read(*,*) alfa,beta
c      alfa=1.d0
c      beta=5.d0
      write(*,*) 'zadaj uroven sumu'
      read(*,*) ss
      if (alfa.EQ.beta) then
         write(*,*) 'tento integral sa rovna nule'
         stop
      end if
      if (alfa.GT.beta) then
         pom=alfa
         alfa=beta
         beta=pom
      end if
      if (alfa.EQ.0.d0) alfa=alfa+1.d-3
      step=(beta-alfa)/(dble(nb)-1.d0)
      pom=alfa
      do i=1,11
         x(i)=pom
         xs(i)=real(x(i))
         rr=rand()
         y(i)=a*exp(-b*pom)+ss*(2.0d0*rr-1.0d0)
         ys(i)=real(y(i))
         z(i)=log(y(i))
         sigma(i)=0.0d0
         pom=pom+step
      end do
      call polfit(x,z,sigma,nb,2,0,aa,chi) ! expon. fit
      write(*,*) nb,exp(aa(1)),aa(2)
      call polfit(x,y,sigma,nb,3,0,bb,chi) ! kvadrat. fit
      a11=exp(aa(1))
      b1=-aa(2)   
      ai0=-a*exp(-b*beta)/b+a*exp(-b*alfa)/b
      ai1=-a11*exp(-b1*beta)/b1+a11*exp(-b1*alfa)/b1
      a3=bb(3)
      a2=bb(2)
      a1=bb(1)
      beta3=beta*beta*beta
      beta2=beta*beta
      beta1=beta
      alfa3=alfa*alfa*alfa
      alfa2=alfa*alfa
      alfa1=alfa
      ai2=a3*(beta3-alfa3)/3.d0+a2*(beta2-alfa2)/2.d0+a1*(beta-alfa)
      call metafl('cons')
      call disini()
      call height(24)
c      call hname(50)
      call labdig(2,'X')
      call labdig(2,'Y')
      call name('x-ova os','X')
      call name('y-ova os','Y')
      ystep=0.1*(ys(1)-ys(10))
      xstep=0.1*real(beta-alfa)
      as=real(alfa)
      bs=real(beta)
      call graf(as,bs,as,xstep,ys(10),1.2*ys(1),ys(10)-1.d0,ystep)
      call linwid(4)
      call color('red')
c      call hsymbl(20)
      do i=1,10
       ixs=nxposn(xs(i))
       iys=nyposn(ys(i))
       call symbol(21,ixs,iys)
      end do
      xstep=0.01*(bs-as)
      ystep=0.01*(ys(1)-ys(10))
      as=xs(1)
      do i=1,100
        xs1(i)=as
        ys1(i)=real(a11)*exp(-real(b1)*xs1(i))
        ys2(i)=real(a3*xs1(i)*xs1(i)+a2*xs1(i)+a1)
        as=as+xstep
c        write(*,*) xs1(i),ys1(i),ys2(i)
      end do
c      stop
      call linwid(4)
      call color('blue')
      call curve(xs1,ys1,100)
c      read(*,*) dummy
      call color('yellow')
      call curve(xs1,ys2,100)
      call disfin
c      call curve(xx,yy,n)
      write(*,*) 'hodnota integralu'
      write(*,*) ' presna          z exponencialy        z paraboly'
      write(*,20) ai0,ai1,ai2
 20   format(2X,3(F10.7,9X))
      end