pro read_snap_all, files, num, folder, Part

;+
;   PURPOSE:
;     Use this to read output files from GADGET
; 
;   INPUTS:
;     files	Number of files per snapshot
;     num		Number of snapshot
;     folder	Folder containing the output files
;     
;   OUTPUT:
;     Part	Struct containing the information of all particles
;   
;   AUTHOR:
;     Felipe Garrido Goicovic
;     Instituto de Astrofísica, PUC, Chile
;     fgarrido@astro.puc.cl
; 
;   September, 2013
  

base=folder

FILES=files


Ntot=0L
NtotGas=0L
Ntotwitmass=0L

for fi=0,FILES-1 do begin	

    exts='000'
    exts=exts+strcompress(string(num),/remove_all)
    exts=strmid(exts,strlen(exts)-3,3)
    
;   if FILES gt 1 then begin 
;     f=base+"snapshot_"+exts+"."+string(fi)
;   endif else begin
;     f=base+"snapshot_"+exts
;   endelse
;     f=strcompress(f,/remove_all)
;     print,"reading ",f,"..."
    
    if FILES gt 1 then begin 
      f=base+strcompress("snapshot_"+exts+"."+string(fi),/remove_all)
    endif else begin
      f=base+strcompress("snapshot_"+exts,/remove_all)
    endelse
    print,"reading ",f,"..."

    npart=lonarr(6)	
    massarr=dblarr(6)
    time=0.0D
    redshift=50.0D
    bytesleft=256-6*4 - 6*8 - 8 - 8
    
    la=intarr(bytesleft/2)

    openr,1,f,/f77_unformatted

    readu,1,npart,massarr,time,redshift,la

    N=total(npart)

    pos_t=fltarr(3,N)
    vel_t=fltarr(3,N)
    id_t=lonarr(N)
    mass_t=fltarr(N)
    type_t=intarr(N)
    
    acc_t=fltarr(3,N)
    gacc_t=fltarr(3,N)
    hacc_t=fltarr(3,N)

    dt_t=fltarr(N)

    ind=where((npart gt 0) and (massarr eq 0)) 
    if ind(0) ne -1 then begin
	Nwithmass= total(npart(ind))
        m_t=fltarr(Nwithmass)
    endif else begin	
        Nwithmass= 0
    endelse

    readu,1,pos_t
    readu,1,vel_t
    readu,1,id_t

    if Nwithmass gt 0 then begin
      readu,1,m_t
    endif

    NGas=npart(0)
    NHalo=npart(1)
    NDisk=npart(2)
    NBulge=npart(3)
    NStars=npart(4)
    NBound=npart(5)

    if Ngas gt 0 then begin
        u_t=fltarr(Ngas)
        readu,1,u_t

        rho_t=fltarr(Ngas)
        readu,1,rho_t
        
        h_t=fltarr(Ngas)
        readu,1,h_t
    endif
    
    if NBound gt 0 then begin
        bh_data=fltarr(NBound)
        readu,1,bh_data
        readu,1,bh_data
        readu,1,bh_data
    endif
    
    readu,1,acc_t
    readu,1,gacc_t
    readu,1,hacc_t

    readu,1,dt_t

    close,1
    
    if Ngas gt 0 then type_t(0:Ngas-1)=0
    if Nhalo gt 0 then type_t(Ngas:Ngas+Nhalo-1)=1
    if Ndisk gt 0 then type_t(Ngas+Nhalo:Ngas+Nhalo+Ndisk-1)=2
    if Nbulge gt 0 then type_t(Ngas+Nhalo+Ndisk:Ngas+Nhalo+Ndisk+Nbulge-1)=3
    if Nstars gt 0 then type_t(Ngas+Nhalo+Ndisk+Nbulge:Ngas+Nhalo+Ndisk+Nbulge+Nstars-1)=4
    if Nbound gt 0 then type_t(Ngas+Nhalo+Ndisk+Nbulge+Nstars:Ngas+Nhalo+Ndisk+Nbulge+Nstars+Nbound-1)=5
    
    NGas_wm=0.0D
    NHalo_wm=0.0D
    NDisk_wm=0.0D
    NBulge_wm=0.0D
    NStars_wm=0.0D
    NBound_wm=0.0D
    
    
    if Ngas gt 0 then begin
      if massarr(0) eq 0 then begin
        Ngas_wm=Ngas
	mass_t(0:Ngas-1)=m_t(0:Ngas-1)
      endif else begin
	mass_t(0:Ngas-1)=massarr(0)
      endelse
    endif
    if Nhalo gt 0 then begin
      if massarr(1) eq 0 then begin
        Nhalo_wm=Nhalo
	mass_t(Ngas:Ngas+Nhalo-1)=m_t(Ngas_wm:Ngas_wm+Nhalo-1)
      endif else begin
	mass_t(Ngas:Ngas+Nhalo-1)=massarr(1)
      endelse
    endif
    if Ndisk gt 0 then begin
      if massarr(2) eq 0 then begin
        Ndisk_wm=Ndisk
	mass_t(Ngas+Nhalo:Ngas+Nhalo+Ndisk-1)=m_t(Ngas_wm+Nhalo_wm:Ngas_wm+Nhalo_wm+Ndisk-1)
      endif else begin
	mass_t(Ngas+Nhalo:Ngas+Nhalo+Ndisk-1)=massarr(2)
      endelse
    endif
    if Nbulge gt 0 then begin
      if massarr(3) eq 0 then begin
        Nbulge_wm=Nbulge
	mass_t(Ngas+Nhalo+Ndisk:Ngas+Nhalo+Ndisk+Nbulge-1)=m_t(Ngas_wm+Nhalo_wm+Ndisk_wm:Ngas_wm+Nhalo_wm+Ndisk_wm+Nbulge-1)
      endif else begin
	mass_t(Ngas+Nhalo+Ndisk:Ngas+Nhalo+Ndisk+Nbulge-1)=massarr(3)
      endelse
    endif
    if Nstars gt 0 then begin
      if massarr(4) eq 0 then begin
        Nstars_wm=Nstars
	mass_t(Ngas+Nhalo+Ndisk+Nbulge:Ngas+Nhalo+Ndisk+Nbulge+Nstars-1)=$
	m_t(Ngas_wm+Nhalo_wm+Ndisk_wm+Nbulge_wm:Ngas_wm+Nhalo_wm+Ndisk_wm+Nbulge_wm+Nstars-1)
      endif else begin
	mass_t(Ngas+Nhalo+Ndisk+Nbulge:Ngas+Nhalo+Ndisk+Nbulge+Nstars-1)=massarr(4)
      endelse
    endif
    if Nbound gt 0 then begin
      if massarr(5) eq 0 then begin
        Nbound_wm=Nbound
	mass_t(Ngas+Nhalo+Ndisk+Nbulge+Nstars:Ngas+Nhalo+Ndisk+Nbulge+Nstars+Nbound-1)=$
	m_t(Ngas_wm+Nhalo_wm+Ndisk_wm+Nbulge_wm+Nstars_wm:Ngas_wm+Nhalo_wm+Ndisk_wm+Nbulge_wm+Nstars_wm+Nbound-1)
      endif else begin
	mass_t(Ngas+Nhalo+Ndisk+Nbulge+Nstars:Ngas+Nhalo+Ndisk+Nbulge+Nstars+Nbound-1)=massarr(5)
      endelse
    endif
    
    Ntot=Ntot+N
    NtotGas=NtotGas+Ngas
    Ntotwitmass=Ntotwitmass+Nwithmass
    
    if fi eq 0 then begin 
      pos=pos_t
      vel=vel_t
      id=id_t
      mass=mass_t
      type=type_t
      if Ngas gt 0 then begin
	u=u_t
	rho=rho_t
	h=h_t
      endif
      acc=acc_t
      gacc=gacc_t
      hacc=hacc_t
      dt=dt_t
    endif else begin
      pos=[[pos],[pos_t]]
      vel=[[vel],[vel_t]]
      id=[id,id_t]
      mass=[mass,mass_t]
      type=[type,type_t]
      if Ngas gt 0 then begin
	u=[u,u_t]
	rho=[rho,rho_t]
	h=[h,h_t]
      endif
      acc=[[acc],[acc_t]]
      gacc=[[gacc],[gacc_t]]
      hacc=[[hacc],[hacc_t]]
      dt=[dt,dt_t]
    endelse
  
endfor

Time_snap=replicate(time,Ntot)

if NtotGas eq 0 then begin
  Part=CREATE_STRUCT('Pos',fltarr(3,Ntot),'Vel',fltarr(3,Ntot),'Id',lonarr(Ntot),'Mass',fltarr(Ntot),'Type',intarr(Ntot),'Acc',fltarr(3,Ntot),'Gacc',fltarr(3,Ntot),'Hacc',fltarr(3,Ntot),'Dt',dblarr(Ntot),'Time',fltarr(Ntot))
endif else begin
  Part=CREATE_STRUCT('Pos',fltarr(3,Ntot),'Vel',fltarr(3,Ntot),'Id',lonarr(Ntot),'Mass',fltarr(Ntot),'Type',intarr(Ntot),'U',dblarr(NtotGas),'Rho',dblarr(NtotGas),'Acc',fltarr(3,Ntot),'Gacc',fltarr(3,Ntot),'Hacc',fltarr(3,Ntot),'Dt',dblarr(Ntot),'Time',fltarr(Ntot)) 
endelse

Part.Pos=pos
Part.Vel=vel
Part.Id=id
Part.Mass=mass
Part.Type=type

if NtotGas ne 0 then begin
Part.U=u
Part.Rho=rho
endif

Part.Acc=acc
Part.Gacc=gacc
Part.Hacc=hacc

Part.Dt=dt

Part.Time=Time_snap

end