Python FortranFile.read_ints方法代码示例

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def read_output(path, header_only=True):
    f = FortranFile(path, 'r')
    ncpu = f.read_ints()
    dim = f.read_ints()
    nparts = f.read_ints()
    if header_only:
        f.close()
        return ncpu, dim, nparts
    f.read_ints()
    f.read_ints()
    f.read_ints()
    f.read_ints()
    f.read_ints()

    x = f.read_reals(dtype=np.float64)
    y = f.read_reals(dtype=np.float64)
    z = f.read_reals(dtype=np.float64)

    vx = f.read_reals(dtype=np.float64)
    vy = f.read_reals(dtype=np.float64)
    vz = f.read_reals(dtype=np.float64)

    m = f.read_reals(dtype=np.float64)

    part_ids = f.read_ints()

    birth = f.read_reals(dtype=np.float32)

    f.close()
    return ncpu, dim, nparts, x, y, z, part_ids

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def read_association(listfile):
    assocFile = FortranFile(listfile, 'r')
    nassoc, columns = assocFile.read_ints()
    _tmp = (assocFile.read_reals(dtype=np.float32)).reshape((columns, nassoc)).transpose()
    assoc = pd.DataFrame(_tmp,
                         columns=['halo_id', 'level', 'halo_mass', 'gal_id', 'gal_mass'])

    assoc[['halo_id', 'level', 'gal_id']] =  assoc[['halo_id', 'level', 'gal_id']].astype(np.int32)
    return assoc

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def read_halo_list(listfile):
    haloFile = FortranFile(listfile, 'r')
    nhalos, columns = haloFile.read_ints()
    _tmp = (haloFile.read_reals(dtype=np.float32)).reshape((columns, nhalos)).transpose()
    halos = pd.DataFrame(_tmp,
                         columns=['id', 'level', 'mass', 'x', 'y', 'z', 'rvir'])
    halos[['id', 'level']] = halos[['id', 'level']].astype(int)

    return halos

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
    def tmp():
        ff = FortranFile(filename)
        h = {}
        h["nbodies"] = ff.read_ints()
        h["massp"] = ff.read_ints()
        h["aexp"] = ff.read_reals(dtype=np.int32)
        h["omega_t"] = ff.read_reals(dtype=np.int32)
        h["age_univ"] = ff.read_reals(dtype=np.int32)
        h["n_halos"], h["n_subhalos"] = ff.read_ints()

        for i in tqdm(range(h["n_halos"] + h["n_subhalos"])):
            infos = {
                "header": h
            }
            infos["nparts"] = ff.read_ints()
            infos["members"] = ff.read_ints()
            infos["idh"] = ff.read_ints()
            infos["timestep"] = ff.read_ints()
            infos["hlevel"], infos["hosthalo"], infos["hostsub"], infos["nbsub"], infos["nextsub"] = ff.read_ints()

            infos["mhalo"] = ff.read_reals(dtype=np.int32)
            infos["pos"] = ff.read_reals(dtype=np.int32)
            infos["speed"] = ff.read_reals(dtype=np.int32)
            infos["L"] = ff.read_reals(dtype=np.int32)
            infos["r"], infos["a"], infos["b"], infos["c"] = ff.read_reals(dtype=np.int32)
            infos["ek"], infos["ep"], infos["et"] = ff.read_reals(dtype=np.int32)
            infos["spin"] = ff.read_reals(dtype=np.int32)
            if not dm_type:
                ff.read_reals()
            infos["rvir"],infos["mvir"], infos["tvir"], infos["cvel"] = ff.read_reals(dtype=np.int32)
            ff.read_reals()
            if not dm_type:
                infos["npoints"] = ff.read_ints()
                infos["rdum"] = ff.read_reals(dtype=np.int32)
                infos["density"] = ff.read_reals(dtype=np.int32)

            if low_mem != None:
                try:
                    keys = list(low_mem)
                except:
                    keys = ['nparts', 'members']

                tmp = {}
                for key in keys:
                    try:
                        tmp[key] = infos[key]
                    except KeyError:
                        print('Invalid key {}, can be any of', infos['keys'])
                yield tmp
            else:
                yield infos
        ff.close()

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
 def read_fortran_FFTfield(self):
       """
       Read a fortran binary file from FFTW assert all Nyquist
       entries to be real.
       """
       f=FortranFile(self.infile,'r')
       Ng=f.read_ints(dtype=np.int32)[0]
       print('Fortran file Ngrid='+str(Ng))
       if (Ng != self.Ngrid):
             print('Ngrid values are not equal!')
       dcr=f.read_reals(dtype=np.complex64)
       dcr=np.reshape(dcr,(Ng//2+1,Ng,Ng),order='F')
       return dcr

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def read_fortran_FFTfield(infile):
    """
    Read a Half-Field with FFTW indexing from
    a Fortran Unformatted Binary file. The first
    entry is a single integer.
    """
    f=FortranFile(infile,'r')
    Ngrid=f.read_ints(dtype=np.int32)[0]
    print('Fortran file Ngrid='+str(Ngrid))
    dcr=f.read_reals(dtype=np.complex64)
    dcr=np.reshape(dcr,(Ngrid//2+1,Ngrid,Ngrid),order='F')
    dcr.dump(infile+'.pickle') # Save infile as a pickle
    return dcr

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def particles_in_halo(tree_brick, start=0, end=None, fun_filter=lambda x: True):
    ''' Open a tree bricks file and associate to each halo the corresponding particles.
    '''
    # Open file
    f = FortranFile(tree_brick, 'r')

    # Give a value to end, by default start + 1
    if end == None:
        end = start + 1

    # Read headers
    nbodies = f.read_ints()[0]
    f.read_reals(dtype=np.float32)
    aexp = f.read_reals(dtype=np.float32)
    f.read_reals(dtype=np.float32)
    age = f.read_reals(dtype=np.float32)
    nhalo, nsubhalo = f.read_ints()
    halo_tot = nhalo + nsubhalo

    halos = {}
    for i in tqdm(range(halo_tot)):
        parts = f.read_ints()[0]
        members = f.read_ints()
        this_id = f.read_ints()[0]
        if (start <= this_id and this_id < end and fun_filter(this_id)):
            for dm_particle_id in members:
                if not halos.has_key(this_id):
                    halos[this_id] = []

                halos[this_id].append(dm_particle_id)
        elif this_id >= end:
            break
        f.read_ints()

        # Irrelevant
        level, hosthalo, hostsub, nbsub, nextsub = f.read_ints()
        mstar = 1e11 * f.read_reals(dtype=np.float32)
        px, py, pz = f.read_reals(dtype=np.float32)
        f.read_reals(dtype=np.float32)
        f.read_reals(dtype=np.float32)
        rad = f.read_reals(dtype=np.float32)[0]
        f.read_reals(dtype=np.float32)
        f.read_reals(dtype=np.float32)
        rvir, mvir, tvir, cvel = f.read_reals(dtype=np.float32)
        f.read_reals(dtype=np.float32)

    f.close()
    return halos

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def read_galaxy_list(listfile):
    galFile = FortranFile(listfile, 'r')
    print(listfile)
    ngal, columns = galFile.read_ints()
    _tmp = (galFile.read_reals(dtype=np.float32)).reshape((columns, ngal)).transpose()
    galaxies = pd.DataFrame(_tmp,
                            columns=['id', 'vt', 'dvz', 'dvr', 'dvtheta', 'mass', 'x', 'y', 'z'])
    galaxies.id.astype(int)

    galaxies['sigma'] = 1/3.*np.sqrt(galaxies.dvz**2 + galaxies.dvtheta**2 + galaxies.dvr**2)
    galaxies['sigmaoverv'] = galaxies.sigma / galaxies.vt
    galaxies['elliptic'] = galaxies.sigmaoverv > 1.5
    galaxies['spiral'] = galaxies.sigmaoverv < 0.8

    return galaxies

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def read_atomic_data(elements=['H', 'He', 'C',     # twelve most abundant elements
                               'N', 'O', 'Ne',
                               'Mg', 'Si', 'S', 
                               'Ar', 'Ca', 'Fe', ] , 
                     data_directory= 'sunnei/AtomicData',   # not robust!  Works when calling from the directory that sunnei is in
                     screen_output=False):

    '''
    This routine reads in the atomic data to be used for the
    non-equilibrium ionization calculations.
 
    Instructions for generating atomic data files
    =============================================
    
    The atomic data files are generated from the routines described by
    Shen et al. (2015) and are available at:
    
    https://github.com/ionizationcalc/time_dependent_fortran
    
    First, run the IDL routine 'pro_write_ionizrecomb_rate.pro' in the
    subdirectory sswidl_read_chianti with optional parameters: nte
    (number of temperature bins, default=501), te_low (low log
    temperature, default=4.0), and te_high (high log temperature,
    default=9.0) to get an ionization rate table.  The routine outputs
    the file "ionrecomb_rate.dat" which is a text file containing the
    ionization and recombination rates as a function of temperature.
    This routine requires the atomic database Chianti to be installed
    in IDL.

    Second, compile the Fortran routine 'create_eigenvmatrix.f90'.
    With the Intel mkl libraries it is compiled as: "ifort -mkl
    create_eigenvmatrix.f90 -o create.out" which can then be run with
    the command "./create.out".  This routine outputs all the
    eigenvalue tables for the first 28 elements (H to Ni).

    As of 2016 April 7, data from Chianti 8 is included in the
    CMEheat/AtomicData subdirectory.
    '''

    if screen_output:
        print('read_atomic_data: beginning program')
    
    from scipy.io import FortranFile

    '''
    Begin a loop to read in the atomic data files needed for the
    non-equilibrium ionization modeling.  The information will be
    stored in the atomic_data dictionary.

    For the first element in the loop, the information that should be
    the same for each element will be stored at the top level of the
    dictionary.  This includes the temperature grid, the number of
    temperatures, and the number of elements.

    For all elements, read in and store the arrays containing the
    equilibrium state, the eigenvalues, the eigenvectors, and the
    eigenvector inverses.
    '''

    atomic_data = {}
    
    first_element_in_loop = True

    for element in elements:

        if screen_output:
            print('read_atomic_data: '+element)

        AtomicNumber = AtomicNumbers[element]
        nstates = AtomicNumber + 1

        filename = data_directory + '/' + element.lower() + 'eigen.dat'
        H = FortranFile(filename, 'r')

        nte, nelems = H.read_ints(np.int32)
        temperatures = H.read_reals(np.float64)
        equistate = H.read_reals(np.float64).reshape((nte,nstates))
        eigenvalues = H.read_reals(np.float64).reshape((nte,nstates))
        eigenvector = H.read_reals(np.float64).reshape((nte,nstates,nstates))
        eigenvector_inv = H.read_reals(np.float64).reshape((nte,nstates,nstates))
        c_rate = H.read_reals(np.float64).reshape((nte,nstates))
        r_rate = H.read_reals(np.float64).reshape((nte,nstates))      
        
        if first_element_in_loop:
            atomic_data['nte'] = nte
            atomic_data['nelems'] = nelems  # Probably not used but store anyway
            atomic_data['temperatures'] = temperatures
            first_element_in_loop = False
        else: 
            assert nte == atomic_data['nte'], 'Atomic data files have different number of temperature levels: '+element
            assert nelems == atomic_data['nelems'], 'Atomic data files have different number of elements: '+element
            assert np.allclose(atomic_data['temperatures'],temperatures), 'Atomic data files have different temperature bins'

        atomic_data[element] = {'element':element,
                                'AtomicNumber':AtomicNumber,
                                'nstates':nstates,
                                'equistate':equistate,
                                'eigenvalues':eigenvalues,
                                'eigenvector':eigenvector,
                                'eigenvector_inv':eigenvector_inv,
#.........这里部分代码省略.........

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
parser.add_argument('--in', dest='infile', type=str, help='Path to the output file of compute_extrema (in hdf format)', required=True)
parser.add_argument('--out', '-o', type=str, help='Prefix of the outputs')
parser.add_argument('--infofile', type=str, help='Path to the information file (the one containing units of RAMSES, …)')

args = parser.parse_args()

# read the info file
infos = dict()
infos['headers'] = pd.read_csv(args.infofile, sep=' *= *', nrows=19, names=['key', 'value'], index_col='key').T
infos['domain'] = pd.read_csv(args.infofile, delim_whitespace=True, skiprows=20)


# read the center
from scipy.io import FortranFile
ff = FortranFile('data/halo_536-centers.bin')
ff.read_ints() # don't care
outputs = ff.read_ints()
centers = ff.read_reals().reshape(len(outputs), 3)
mins    = ff.read_reals().reshape(len(outputs), 3)
span    = ff.read_reals().reshape(len(outputs), 3)
maxs    = ff.read_reals().reshape(len(outputs), 3)

# create the output dir if required
# if not os.path.isdir(args.out):
#     os.mkdir(args.out)


HDF = pd.HDFStore(args.infile)
# read the data
df    = HDF['extremas']
dens  = HDF['dens']

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
#------------------------------------------------------------------------------
# open file and read
#------------------------------------------------------------------------------
# Parameters
path_eigentb = '/Users/ccai/Works/Project/Ionization_calc/Code_develop/\
ionization_code_reu2016/python_script/chianti_8/'
element = 'o'
file_name = element+'eigen.dat'

# Open file
file_eigentb = path_eigentb + file_name
f = FortranFile(file_eigentb, 'r')

# Read file
[nte, natom]=f.read_ints(dtype=np.int32)
te_arr = f.read_reals(dtype=np.float64)
eqistate = f.read_reals(dtype=np.float64).reshape((natom+1, nte), order='F')
eigenvals = f.read_reals(dtype=np.float64).reshape((natom+1, nte), order='F')
eigenvector = f.read_reals(dtype=np.float64).reshape((natom+1, natom+1, nte), order='F')
eigenvector_invers = f.read_reals(dtype=np.float64).reshape((natom+1, natom+1, nte), order='F')

# Close file
f.close()


# Note from Nick: I copied the next three routines to time_advance.py
# but am leaving these also here for now.


#------------------------------------------------------------------------------

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def main(sourcedir, sourcefile, savedir, saveheader, moving=0, incr=1, imin=0, imax=-1 ):
    """
    imin --> minimum iteration above which files will be written
    incr --> iteration interval to write files at
    """

    # ogdir = os.getcwd()
    # os.chdir(savedir)

    MakeOutputDir(savedir)

    #link in source file
    if not os.path.isfile('{}/{}'.format(savedir, sourcefile)):
        # cmd('ln -s {}/{} {}'.format(sourcedir, sourcefile, sourcefile))
        cmd('ln -s {}/{} {}/{}'.format(sourcedir, sourcefile, savedir, sourcefile))


    #GET GRID/FILE DIMENSIONS
    f = FortranFile('{}/{}'.format(savedir, sourcefile), 'r')
    dims = f.read_ints('uint32')[2:6]
    #Number of iterations, Number of points in each direction, number of parameters, needed to read whole file
    nj, nk, nl, nq = dims

    #FORMATS IN BC201 FILE
    bc201formats = [
                        ('ints', '<i4', (1,7)), #IGRID,ISTEP,NJ,NK,NL,NQ,NQC
                        ('tvr', 'float64'),     #TVREF
                        ('dtr', 'float64'),     #DTVREF
                        ('xyz', 'float64', (3,nj,nk,nl)), #XYZ points
                        ('q', 'float64', (nq,nj,nk,nl)),  #Q variables (nq variables, njXnkXnl values for each)
                        ('iblank', 'int32', (nj,nk,nl)),  #IBLANK status of each grid point
                    ]

    #########################
    #Save pressure coeff history at certian locations
    #! I want to write out a subset through time (every iteration!)
    ports = [8,13,30,45]
    out = open(savedir+"/cp_history.dat","w")
    out.write("ITER")
    for p in ports:
        out.write(" {}_x {}_cp".format(p,p))
    out.write("\n")
    #############################


    #RESTART FILE AND READ EACH RECORD
    f = FortranFile('{}/{}'.format(savedir, sourcefile), 'r')
    keep_reading = True
    irecord = -1
    while (keep_reading is True):

        irecord += 1

        #READ CURRENT RECORD
        b = f.read_record(bc201formats)



        #GET CURRENT ITERATION
        istep = b['ints'][0][0][1]

        # print(istep)


        #DONT WRITE UNDESIRED FILES (SKIP)
        #only write files at specified interval
        if (istep % incr != 0):
            continue
        #Don't write files below start iter
        if (istep < imin):
            continue
        #Stop reading when max desired iteration is reached
        if istep > imax - incr:
            keep_reading = False




        #WRITE GRID POINTS TO PLOT3D FORMATTED FILE

        if moving:
            #Grid is moving, write to file at each iteration
            savename = '{}/{}.x.{}'.format(savedir, saveheader, istep)
            WriteGrid(savename, b, nj, nk, nl)
        elif irecord == 0:
           #Grid is stationary, write to file only once
            savename = '{}/{}.x'.format(savedir, saveheader)
            WriteGrid(savename, b, nj, nk, nl)


        #CALCULATE CP FOR EACH POINT
        cps = np.zeros((nj,nk,nl))
        # print(cps.shape)
        for j in range(nj):
            for k in range(nk):
                for l in range(nl):
                    #print(j,k,l,b['q'].shape)
                    q = np.zeros(nq)
                    for iq in range(nq):
                        q[iq] = b['q'][0][iq][j][k][l]
#.........这里部分代码省略.........

# 需要导入模块: from scipy.io import FortranFile [as 别名]
# 或者: from scipy.io.FortranFile import read_ints [as 别名]
def read():
	f=FortranFile(BinaryData,'r')
	n_vals=f.read_record('i4,i4,i4,i4,i4,(12,)i4')
	d_vals=f.read_reals('f4')
	xyz_vals=f.read_record('(500,)f4,(500,)f4,(500,)f4')
	a_vals=f.read_reals('i8')
	ndr_vals=f.read_ints('i4')
	ag_vals=f.read_record('(3,)i8')
	f.close()
	return n_vals,d_vals,xyz_vals,a_vals,ndr_vals,ag_vals