program Polovina_diska

implicit none

real*8:: x,ran2, y,duzina,sumx,sumy ! x i y koordinate diska, duzina rastojanje tacke od koordinatnog pocetka
integer::i,idum,idum1,ibr !i Korak Brojaca
parameter R=0.25,k=11000 !Poluprecnik diska 
	 real*8, parameter::a1=0., b1=0.5,a2=-0.5,b2=0.5,pi=3.142

open(30,file='centar.dat',status='unknown') !Otvaranje datoteke u koju se smestaj fajlovi

idum=-21908
idum1=-21312
sumx=0
sumy=0
duzina=0 !Pocetna suma nasumicnih centara masa
ibr=0
do i=1,k !Brojac koji broji korake

x=ran2(idum) !Nasumicna vrednost broja x
y=(ran2(idum1)-0.5) !Nasumicna vrednost broja y

if (((x*x+y*y).le.0.25)) then !If petlja koja ispituje da li tacka pripada poudisku
ibr=ibr+1
duzina=duzina+(sqrt(x*x+y*y)) !Sumiranje duzina. Rezultat ukupna duzina svih centara masa koji pripadaju disku
sumx=sumx+x*(b1-a1)*(b2-a2)*2
sumy=sumy+y*(b1-a1)*(b2-a2)*2

write (30,100) x,y,sqrt(x*x+y*y),sumx/ibr,sumy/ibr,duzina/ibr ! x,y koordinata, vektorsko rasrojanje od centra, srednja vrednost centra masa

endif
enddo

print*, sumx/ibr,sumy/ibr,duzina/ibr

100 format(10(e13.7,2x))

stop

end




!     *****************************RAN2****************************
      function ran2(idum)
      integer*4 idum,IM1,IM2,INM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      real*8 ran2,AM,EPS,RNMX
      parameter(IM1=2147483563,IM2=2147483399,AM=1.d0/IM1,INM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,NTAB=32,NDIV=1+INM1/NTAB,EPS=1.2d-7,RNMX=1.-EPS)
!     Long period (>2x10^{18}) random number generator of L'Ecuyer
!     with Bays-Durham shuffle and added safeguards. Returns a uniform
!     random deviate between 0.0 and 1.0 (exclusive of the endpoint
!     values). Call with idum a negative integer to initialise;
!     thereafter, do not alter idum between successive deviate in a
!     sequence. RNMX should approximate the largest floating value
!     that is less than 1.
      integer*4 idum2,j,k,iv(NTAB),iy
      save iv,iy,idum2
      data idum2/123456789/,iv/NTAB*0/,iy/0/

      if(idum .le. 0) then             !Initialise
         idum = max(-idum,1)           !Be sure to prevent idum=0
         idum2 = idum
         do 100, j=NTAB+8,1,-1         !Load the shuffle table (after
            k=idum/IQ1                 ! 8 warm-ups)
            idum = IA1*(idum-k*IQ1)-k*IR1
            if(idum .lt. 0) idum = idum + IM1
            if(j .le. NTAB) iv(j) = idum
 100     continue
         iy = iv(1)
      end if
      k = idum/IQ1                     !Start here when not initialising
      idum = IA1*(idum-k*IQ1)-k*IR1    !Compute idum=mod(IA1*idum,IM1)
                                       !without overflow by Schrage's
                                       !method
      if(idum .lt. 0) idum = idum + IM1
      k = idum2/IQ2
      idum2 = IA2*(idum2-k*IQ2)-k*IR2  !Compute idum2=mod(IA2*idum2,IM2)
                                       !likewise
      if(idum2 .lt. 0) idum2 = idum2 + IM2
      j = 1 + iy/NDIV                  !Will be in the range 1:NTAB
      iy = iv(j) - idum2               !Here idum is shuffled, idum and
      iv(j) = idum                     !idum2 are combined to generate
                                       !output
      if(iy .lt. 1) iy = iy + INM1
      ran2 = min(AM*iy,RNMX)           !Because users do not expect
                                       !endpoint values
      return
      end
!     ******************************************************************







