forked from ilkkavir/EFISR
-
Notifications
You must be signed in to change notification settings - Fork 0
/
readParametersFromFileGUISDAP.m
157 lines (131 loc) · 5.45 KB
/
readParametersFromFileGUISDAP.m
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
147
148
149
150
151
152
153
154
155
156
function param = readParametersFromFileGUISDAP( fpath )
%function [Ne,NeStd,Ti,TiStd,Tr,TrStd,Coll,CollStd,Vi,ViStd,Comp,CompStd,status,chisqr,ts,te,mlt,llhT,llhR,azelR,r,h,phi,site,ecefS,llhgS,llhmS,kS,B] = readParametersFromFileGUISDAP( fpath )
%
% [Ne,NeStd,Ti,TiStd,Tr,TrStd,Coll,CollStd,Vi,Vistd,Comp,CompStd,status,
% chisqr,ts,te,mlt,llhT,llhR,azelR,r,h,phi,site,ecefS,llhgS,llhmS,kS,B]
% = readParametersFromFileGUISDAP( fpath )
%
% Read plasma parameters and related information from GUISDAP fit
% result files. This function reads exactly one file
%
% INPUT:
% fpath path to the GUISDAP output file
%
% OUTPUT:
%
% param a struct with fields:
%
% Ne Electron density (m^-3)
% NeStd Ne standard deviation (m^-3)
% Ti Ion temperature (K)
% TiStd Ti standard deviation (K)
% Tr Temperature ratio (electron/ion) (dimensionless)
% TrStd Tr standard deviation
% Coll Ion-neutral collision frequency (s^-1)
% CollStd Coll standard deviation (s^-1)
% Vi Line-of-sight ion velocity (ms^-1)
% ViStd Vi standard deviation (ms^-1)
% status fit status (0 for success)
% chisqr chi-squared
% ts integration start time (unix time)
% te integration end time (unix time)
% llhT latitude, longitude, height of the transmitter (deg,deg,m)
% llhR latitude, longitude, height of the receiver (deg,deg,m)
% azelR azimuth and elevation angles of the receiver (deg,deg)
% r range (m)
% h height (m)
% phi scattering angle (deg)
% site EISCAT site name (char)
%
%
% the following parameters have been removed from the output on 18 Aug 2021
% ecefS location of the scattering volume(s) in cartesian ecef
% (earth-centred earth-fixed) coordinates (m)
% llhgS geodetic coordinates of the scattering volume(s). (deg,deg,m)
% llhmS altitude-adjusted corrected geomagnetic (v2)
% coordinates of the scattering volume (deg,deg,m). The
% last coordinate is geocentric height (radial distance
% - Re)
% kS scattering wave vector direction(s) (unit vector(s)) in
% cartesian ecef coordinates
% B magnetic field vector in local North-East-Down
% coordinates (nT)
% mlt magnetic local time at integration end time (hours)
%
%
% IV 2016
%
% read the data
load(fpath);
param = struct();
% parameters and their standard deviations
param.Ne = r_param(:,1);
param.Ti = r_param(:,2);
param.Tr = r_param(:,3);
param.Coll = r_param(:,4);
param.Vi = r_param(:,5);
param.Comp = r_param(:,6);
param.NeStd = r_error(:,1);
param.TiStd = r_error(:,2);
param.TrStd = r_error(:,3);
param.CollStd = r_error(:,4);
param.ViStd = r_error(:,5);
param.CompStd = r_error(:,6);
% fit status
param.status = r_status;
% chi-squared
param.chisqr = r_res(:,1);
% number of velocity points
nv = length(param.Vi);
% integration start and end time as unix time
param.ts = repmat( posixtime(datetime(r_time(1,1),r_time(1,2),r_time(1,3),r_time(1,4), ...
r_time(1,5),r_time(1,6))) , nv , 1 );
param.te = repmat( posixtime(datetime(r_time(2,1),r_time(2,2),r_time(2,3),r_time(2,4), ...
r_time(2,5),r_time(2,6))) , nv , 1 );
% transmitter location (lat,lon,height)
param.llhT = repmat( r_XMITloc , nv , 1 );
param.llhT(:,3) = param.llhT(:,3)*1000; % height is in km in guisdap files
% receiver location (lat,lon,height)
param.llhR = repmat( r_RECloc , nv , 1 );
param.llhR(:,3) = param.llhR(:,3)*1000; % height is in km in guisdap files
% azimuth and elevation of the receiver antenna
param.azelR = repmat( [ r_az , r_el ] , nv , 1 );
% range (distance from target to receiver? check this!)
param.r = r_range.*1000;
% height as calculated by GUISDAP
param.h = r_h.*1000;
% scattering angle calculated by GUISDAP
param.phi = repmat( r_SCangle.*180./pi , nv , 1);
% site name
param.site = repmat( name_site , nv , 1 );
% % scattering volume locations etc.
% param.ecefS = zeros(nv,3);
% param.llhgS = zeros(nv,3);
% param.llhmS = zeros(nv,3);
% param.kS = zeros(nv,3);
% param.B = zeros(nv,3);
% param.mlt = zeros(nv,1);
% % should calculate only unique values like in Daedalus/calculateCoordinatesMadrigal.m!
% % this is now done for hdf5 data in readVelocitiesGUISDAP. Could remove the following lines and work also with matlab file data there
% for hh = 1:nv
% [param.ecefS(hh,:),param.llhgS(hh,:),param.llhmS(hh,:),param.kS(hh,:)] = ...
% systemGeometry2scatteringGeometry(param.llhT(hh,:),param.llhR(hh,:), ...
% param.azelR(hh,:),param.r(hh), ...
% datetime(param.ts(hh,:), ...
% 'ConvertFrom', ...
% 'posixtime'));
% % magnetic field direction from igrf
% % maximum altitude of igrf is 600 km...
% if param.llhgS(hh,3) <= 600000
% try
% param.B(hh,:) = igrfmagm( param.llhgS(hh,3) , param.llhgS(hh,1) , param.llhgS(hh,2), ...
% decyear(r_time(2,:)));
% catch
% param.B(hh,:) = igrfmagm( param.llhgS(hh,3) , param.llhgS(hh,1) , param.llhgS(hh,2), ...
% 2020);
% end
% % magnetic local time
% param.mlt(hh) = magneticLocalTime(datetime(param.ts(hh,:),'ConvertFrom','posixtime'),param.llhmS(hh,2));
% end
% end
end