C ______________________________ Analysis of tables with profiles
C		                 for barred galaxies simulations
C		                     
c-------------------------------------------------------------------- 
      INCLUDE  'PMgalaxy2.h'
      REAL         INPUT,moment
      Character    FileASCII*100

       write(*,*)' enter 6 digits of the snapshot Example: 0.6235'
       read(*,*) moment

        write(FileASCII,'("RESULTS/Result.SPH.",f6.4,".dat")')
     &    moment
        write (*,*) FileASCII 
      OPEN(37,FILE=FileASCII,STATUS='UNKNOWN') 

         write(FileASCII,'("RESULTS/Result.CYL.",f6.4,".dat")')
     &    moment
      OPEN(39,FILE=FileASCII,STATUS='UNKNOWN') 


      CALL ReadSPH
      CALL FreqSPH


      Return
      End
c-------------------------------------------------------------------- 
c                     Read table SPH
      SUBROUTINE ReadSPH
      INCLUDE 'PMgalaxy2.h'
      Character    Line*100

      Do i=1,2
         Read(37,'(A)') Line
         write(*,*) Line
      EndDo 
      Read(37,*) AEXPN,dA,Time,iStep,Ndisk,Nbulge,aMdisk
      write (*,*)  
      write (*,*)  ' A    =',AEXPN
      write (*,*)  ' Time =',Time,'Gyrs'
      write (*,*)  ' Mdisk=',aMdisk,'Msun'
      Do i=1,5
         Read(37,'(A)') Line
c         write(*,*) Line
      EndDo 

      Ntab =1
 10   read(37,*,err=50,end=50)Radout(Ntab),Radd(Ntab),
     &         Vch(Ntab),Vcdisk(Ntab),Vcbulge(Ntab),
     &         Vctot(Ntab),Nradh(Ntab),Nradhout(Ntab),
     &         Nradd(Ntab),Nradout(Ntab),
     &         Nradb(Ntab),Nradbout(Ntab)
      Ntab =Ntab +1
      goto 10

 50   Ntab =Ntab -1
      write (*,*)  ' Read ',Ntab,' lines from table SPH'

      Return
      End
c-------------------------------------------------------------------- 
c                     Frequencies: Omega,epicycle
      SUBROUTINE FreqSPH
      INCLUDE 'PMgalaxy2.h'

      Do i=2,Ntab-1
         Omega = Vctot(i)/Radout(i)
         dVdRL = (Vctot(i)-Vctot(i-1))/(Radout(i)-Radout(i-1))
         dVdRR = (Vctot(i+1)-Vctot(i))/(Radout(i+1)-Radout(i))
         dVdR  = (dVdRR+dVdRL)/2.
         aKappa= sqrt(2.)*Omega*sqrt(1.+dVdR/Omega)
         write (*,100)i,Radout(i),Vctot(i),Omega,aKappa,Omega-aKappa/2.,
     &      dVdR,
     &      dVdRL,dVdRR
      EndDo 
 100     format(i3,8g12.4)
      Return
      End