%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                       %%
%%  Version JA-4.56 of Feb 2021                                          %%
%%                                                                       %%
%%  This is a matlab routine to load parody G files                      %%
%%  Type fname='your_G_file_name'                                        %%
%%  and then invoke this script                                          %%
%%  output is scaled to dimensional values for the 71% of path model     %%
%%                                                                       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid=fopen(fname,'r');

version=fread(fid,1,'int');
phypar=fread(fid,10,'double')';
time=phypar(1);
DeltaU=phypar(2);
Coriolis=phypar(3);
Lorentz=phypar(4);
Buoyancy=phypar(5);
ForcingU=phypar(6);
DeltaT=phypar(7);
ForcingT=phypar(8);
DeltaB=phypar(9);
ForcingB=phypar(10);

Ek=1/Coriolis;
Ra=Buoyancy*Ek;
Pm=1/DeltaB;
Pr=1/DeltaT;

gridpar=fread(fid,4,'int');
nr=gridpar(1);
nt=gridpar(2);
np=gridpar(3);
azsym=gridpar(4);

r=fread(fid,nr,'double');
theta=fread(fid,nt,'double');
sint=sin(theta);
cost=cos(theta);
phi=[1:np]'/np*2*pi/azsym;

Vr=zeros(np,nt,nr);
Vt=zeros(np,nt,nr);
Vp=zeros(np,nt,nr);
Br=zeros(np,nt,nr);
Bt=zeros(np,nt,nr);
Bp=zeros(np,nt,nr);
T=zeros(np,nt,nr);

a=fread(fid,np*nt*nr,'single');
Vr=reshape(a,np,nt,nr);

a=fread(fid,np*nt*nr,'single');
Vt=reshape(a,np,nt,nr);

a=fread(fid,np*nt*nr,'single');
Vp=reshape(a,np,nt,nr);

a=fread(fid,np*nt*nr,'single');
Br=reshape(a,np,nt,nr);

a=fread(fid,np*nt*nr,'single');
Bt=reshape(a,np,nt,nr);

a=fread(fid,np*nt*nr,'single');
Bp=reshape(a,np,nt,nr);

a=fread(fid,np*nt*nr,'single');
T=reshape(a,np,nt,nr);

fclose(fid);

%scaling to dimensional values

%time in years
ts=300/sqrt(3e-10);
timers=time*ts;

%radius in km
r=r*2260;

% velocity in km/yr
Vr=Vr*2260/ts;
Vt=Vt*2260/ts;
Vp=Vp*2260/ts;

% magnetic field in mT
Br=Br*0.9;
Bt=Bt*0.9;
Bp=Bp*0.9;

% density anomaly field in kg/m^3
% factor Ra*Pm is for expressing dimensionless result in units of \rho\Omega\eta/(g_o D)
% factor 4.25e-8 is Earth value for \rho\Omega\eta/(g_o D)
% constant is arbitrary, adjusted here so that equatorial density anomaly is well
% represented between -2e-4 and 2e-4 kg/m^3

T=T*Ra*Pm*4.25e-8-0.20132;






