-
Notifications
You must be signed in to change notification settings - Fork 37
Expand file tree
/
Copy pathreadgrid.m
More file actions
146 lines (120 loc) · 4.66 KB
/
Copy pathreadgrid.m
File metadata and controls
146 lines (120 loc) · 4.66 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
function xgf=readgrid(inID)
%--------------------------------------------------------
%-----THIS READS A GRID FROM A BINARY FILE CREATED
%-----BY MATLAB OR POSSIBLY FORTRAN (THOUGH THIS IS
%-----NOT YET IMPLEMENTED AS OF 9/15/2016)
%--------------------------------------------------------
narginchk(1,1)
validateattr(inID, {'char'}, {'vector'}, mfilename, 'grid directory', 1)
filename=[inID,filesep,'simsize.dat'];
if ~exist(filename,'file'), error([filename,' does not exist ']), end
fid=fopen(filename,'r');
xgf.lx=fread(fid,3,'integer*4');
fclose(fid);
lx1=xgf.lx(1); lx2=xgf.lx(2); lx3=xgf.lx(3);
lgrid=lx1*lx2*lx3;
lgridghost=(lx1+4)*(lx2+4)*(lx3+4);
gridsize=[lx1,lx2,lx3];
gridsizeghost=[lx1+4,lx2+4,lx3+4];
%%
fin = [inID,'simgrid.dat'];
if ~exist(fin,'file'), error([fin,' does not exist\']), end
fid=fopen(fin,'r');
xgf.x1=fread(fid,lx1+4,'real*8'); %coordinate values
xgf.x1i=fread(fid,lx1+1,'real*8');
xgf.dx1b=fread(fid,lx1+3,'real*8'); %ZZZ - need to check that differences have appropriate ghost cell values, etc.
xgf.dx1h=fread(fid,lx1,'real*8');
xgf.x2=fread(fid,lx2+4,'real*8');
xgf.x2i=fread(fid,lx2+1,'real*8');
xgf.dx2b=fread(fid,lx2+3,'real*8');
xgf.dx2h=fread(fid,lx2,'real*8');
xgf.x3=fread(fid,lx3+4,'real*8');
xgf.x3i=fread(fid,lx3+1,'real*8');
xgf.dx3b=fread(fid,lx3+3,'real*8');
xgf.dx3h=fread(fid,lx3,'real*8');
tmp=fread(fid,lgridghost,'real*8'); %cell-centered metric coefficients
xgf.h1=reshape(tmp,gridsizeghost);
tmp=fread(fid,lgridghost,'real*8');
xgf.h2=reshape(tmp,gridsizeghost);
tmp=fread(fid,lgridghost,'real*8');
xgf.h3=reshape(tmp,gridsizeghost);
tmpsize=[lx1+1,lx2,lx3];
ltmp=prod(tmpsize);
tmp=fread(fid,ltmp,'real*8'); %interface metric coefficients
xgf.h1x1i=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.h2x1i=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.h3x1i=reshape(tmp,tmpsize);
tmpsize=[lx1,lx2+1,lx3];
ltmp=prod(tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.h1x2i=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.h2x2i=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.h3x2i=reshape(tmp,tmpsize);
tmpsize=[lx1,lx2,lx3+1];
ltmp=prod(tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.h1x3i=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.h2x3i=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.h3x3i=reshape(tmp,tmpsize);
%gravity, geographic coordinates, magnetic field strength? unit vectors?
tmp=fread(fid,lgrid,'real*8'); %gravitational field components
xgf.gx1=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8');
xgf.gx2=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8');
xgf.gx3=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8'); %geographic coordinates
xgf.alt=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8');
xgf.glat=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8');
xgf.glon=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8'); %magnetic field strength
xgf.Bmag=reshape(tmp,gridsize);
tmp=fread(fid,lx2*lx3,'real*8'); %magnetic field inclination - only one value for each field line
xgf.I=reshape(tmp,[lx2,lx3]);
tmp=fread(fid,lgrid,'real*8'); %points not to be solved
xgf.nullpts=reshape(tmp,gridsize);
%STUFF PAST THIS POINT ISN'T USED IN FORTRAN CODE BUT INCLUDED IN THE
%GRID FILE FOR COMPLETENESS
if (~feof(fid))
tmpsize=[lx1,lx2,lx3,3];
ltmp=prod(tmpsize);
tmp=fread(fid,ltmp,'real*8'); %4D unit vectors (in cartesian components)
if (feof(fid)) %for whatever reason, we sometimes don't hit eof until after first unit vector read...
return; %lazy as hell
end
xgf.e1=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.e2=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.e3=reshape(tmp,tmpsize);
tmpsize=[lx1,lx2,lx3,3];
ltmp=prod(tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.er=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.etheta=reshape(tmp,tmpsize);
tmp=fread(fid,ltmp,'real*8');
xgf.ephi=reshape(tmp,tmpsize);
tmp=fread(fid,lgrid,'real*8'); %spherical coordinates
xgf.r=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8');
xgf.theta=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8');
xgf.phi=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8'); %cartesian coordinates
xgf.x=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8');
xgf.y=reshape(tmp,gridsize);
tmp=fread(fid,lgrid,'real*8');
xgf.z=reshape(tmp,gridsize);
end
fclose(fid);
end