Bonjour à tous,
J'ai un petit souci avec l'utilisation de parfor.
Quand je fais tourner le code ci-dessous je n'ai aucun problème :
Le calcul s'effectue bien sur deux cœurs et à chaque fois que je passe dans les boucle n et m je calcul bien les variables : lat , lon , SVF , SBVF et SBWF.
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223 tic time_zone = 4 ; %Time zone (i.e. GMT+) ; dst = 0 ; %Daylight saving time ; Lmax = 55.8272 ; %Goegraphic area for calculation ; Lmin = 55.2168 ; lmax = -20.8722 ; lmin = -21.3891 ; step = 0.1 ; LL = Lmin:step:Lmax ; ll = lmin:step:lmax ; R_az = 45 ; %Azimuth angle resolution [°] ; R_di = 0.005 ; %Resolution a direction defined by any azimuth angle [°] ; R_az = R_az * pi / 180 ; %Conversion degre to radian ; matlabpool('open',2); parfor n = 1:length(LL) for m = 1:length(ll) lat = ll(m) ; %Geopoint M(lat,long) ; lon = LL(n) ; display(lat); display(lon); display( (now - time1)*86400 ) ; altM = interp2( lat_long_alt_Reu_STRM(:,:,2),... %Altitude of M [m]; lat_long_alt_Reu_STRM(:,:,1),... lat_long_alt_Reu_STRM(:,:,3),... lon , lat ,'cubic') ; if isnan( altM ) altM = 0 ; end azi = pi:-R_az:-pi ; Hz = zeros(1,length(azi)); for i = 1:length(azi) %Azimuth step angle vector ; k = 1 ; r = 0 ; X = lon ; Y = lat ; while lmin <= Y && Y <= lmax && Lmin <= X && X <= Lmax ; r = r + R_di ; %Geopoint N(X,Y) ; X = lon + r * cos(azi(i)+pi/2) ; Y = lat + r * sin(azi(i)+pi/2) ; altN = interp2( lat_long_alt_Reu_STRM(:,:,2),... %Altitude of N [m] ; lat_long_alt_Reu_STRM(:,:,1),... lat_long_alt_Reu_STRM(:,:,3),... X , Y ,'cubic') ; if isnan( altN ) altN = 0 ; end altMN = altN - altM ; %Altitude difference between M and N [m] ; distMN = vdist( lat , lon , Y , X ) ; %Distance between M and N [m] ; Ele(n,k) = atan( altMN / distMN ) * 180 / pi ; %Elevation of N from M ; k = k + 1 ; end Hz(i) = max(Ele(n,:)) ; %Overshagin of M ; Ele = zeros(1) ; end azi(1) = [] ; Hz(1) = [] ; %Sky View Factor calculation: TotalSky = 0 ; SkyView = 0 ; for i = 1:length(Hz) for j = 1:0.1:90 ; TotalSky = TotalSky + 1 ; if j > Hz(i) SkyView = SkyView + 1 ; end end end SVF = SkyView / TotalSky ; %Sun Beam (and weighted) View Factor calculation: TotalSunPath = 0 ; SunPathView = 0 ; TotalBeamEnergy = 0 ; BeamEnergyView = 0 ; for j=1:365 for i=1:1440 %Note: i represents the minute of the day and j represents the %day of the year. %The equation of time (E): B = ( j - 1 ) * 2 * pi / 365; E = 229.18 * ( 7.5e-5 + 1.868e-3 * cos(B) - 3.2077e-2 * sin(B)- 1.4615e-2 * cos(2*B) - 4.089e-2 * sin(2*B) ); %Convert to west longitude: if lon>0 longitude_standard=360-(time_zone+dst)*15; longitude_in_west_degree=360-lon; else longitude_standard= -(time_zone+dst)*15; longitude_in_west_degree= -lon; end %Difference in minutes between solar time and local time solartime_minus_localtime = 4 *( longitude_standard - longitude_in_west_degree )+ E; %Solar time (minute): solar_time= i + solartime_minus_localtime; if solar_time > 1440 solar_time = solar_time - 1440 ; end if solar_time <= 0 solar_time = 1440 - solar_time ; end %Decimal solar hour solar_hour = solar_time/60 ; %Solar angle calculation: %Hour angle (omega): omega = (solar_hour - 12) * 15 * pi / 180 ; %Declination angle (delta): delta = 23.45 * pi/180 * sin (2*pi*(j+284)/365); %Solar elevation angle (alpha_s): alpha_s = asin( cos(omega) * cos(delta) * cos(lat*pi/180) + sin(delta) * sin(lat*pi/180) ); alpha_s_degree = 180 * alpha_s/pi; % Solar azimuth angle (gamma_s): if lat >= 0 gamma_s = atan2( sin(omega) , cos(omega)*sin(lat*pi/180)-tan(delta)*cos(lat*pi/180) ) + pi ; else gamma_s = atan2( sin(omega) , cos(omega)*sin(lat*pi/180)-tan(delta)*cos(lat*pi/180) ) ; if gamma_s > 0 gamma_s = gamma_s - pi ; else gamma_s = gamma_s + pi ; end end gamma_s_degree = gamma_s * 180 / pi ; if alpha_s_degree > 90 theta_s = ( alpha_s_degree - 90 ) * pi / 180 ; else theta_s = ( 90 - alpha_s_degree ) * pi / 180 ; end if alpha_s > 0 RavR2 = 1.0001 + 0.034221 * cos(B) + 0.001280 * sin(B) + 0.000719 * cos(2*B) + 0.000077 * sin(2*B); %Note: RavR2 is the ratio beteween the mean sun-earth distance %(Rav) and the the actual sun-earth distance depending on the day of the year I0 = 1367 * RavR2 ; %Extraterrestrial radiation in W.m-2 AM = 1 / cos(theta_s) ; BNI = I0*0.7^(AM^0.678) ; TotalSunPath = TotalSunPath + 1 ; TotalBeamEnergy = TotalBeamEnergy + BNI ; ind = find( azi <= gamma_s ) ; if alpha_s_degree > Hz(ind(1)) SunPathView = SunPathView + 1 ; BeamEnergyView = BeamEnergyView + BNI ; end end end end SBVF = SunPathView / TotalSunPath ; SBWVF = BeamEnergyView / TotalBeamEnergy ; end end matlabpool('close'); toc
Maintenant je souhaite récupérer ces variables pour chaque couple (n,m) calculé. Pour cela j'introduis le tableau ViewF, mais Matlab me renvoi l'erreur suivante :
Voici plus bas comment j'utilise ViewF:Error: The variable ViewF in a parfor cannot be classified.
Je ne comprends pas ce qui pose problème.
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237 %*************************************************************************% tic time1 = now ; disp('processing...') time_zone = 4 ; %Time zone (i.e. GMT+) ; dst = 0 ; %Daylight saving time ; Lmax = 55.8272 ; %Goegraphic area for calculation ; Lmin = 55.2168 ; lmax = -20.8722 ; lmin = -21.3891 ; step = 0.1 ; LL = Lmin:step:Lmax ; ll = lmin:step:lmax ; R_az = 45 ; %Azimuth angle resolution [°] ; R_di = 0.005 ; %Resolution a direction defined by any azimuth angle [°] ; R_az = R_az * pi / 180 ; %Conversion degre to radian ; spmd ViewF = zeros(length(LL), length(ll) ,5 ); end matlabpool('open',2); parfor n = 1:length(LL) for m = 1:length(ll) lat = ll(m) ; %Geopoint M(lat,long) ; lon = LL(n) ; display(lat); display(lon); display( (now - time1)*86400 ) ; altM = interp2( lat_long_alt_Reu_STRM(:,:,2),... %Altitude of M [m]; lat_long_alt_Reu_STRM(:,:,1),... lat_long_alt_Reu_STRM(:,:,3),... lon , lat ,'cubic') ; if isnan( altM ) altM = 0 ; end azi = pi:-R_az:-pi ; Hz = zeros(1,length(azi)); for i = 1:length(azi) %Azimuth step angle vector ; k = 1 ; r = 0 ; X = lon ; Y = lat ; while lmin <= Y && Y <= lmax && Lmin <= X && X <= Lmax ; r = r + R_di ; %Geopoint N(X,Y) ; X = lon + r * cos(azi(i)+pi/2) ; Y = lat + r * sin(azi(i)+pi/2) ; altN = interp2( lat_long_alt_Reu_STRM(:,:,2),... %Altitude of N [m] ; lat_long_alt_Reu_STRM(:,:,1),... lat_long_alt_Reu_STRM(:,:,3),... X , Y ,'cubic') ; if isnan( altN ) altN = 0 ; end altMN = altN - altM ; %Altitude difference between M and N [m] ; distMN = vdist( lat , lon , Y , X ) ; %Distance between M and N [m] ; Ele(n,k) = atan( altMN / distMN ) * 180 / pi ; %Elevation of N from M ; k = k + 1 ; end Hz(i) = max(Ele(n,:)) ; %Overshagin of M ; Ele = zeros(1) ; end azi(1) = [] ; Hz(1) = [] ; %Sky View Factor calculation: TotalSky = 0 ; SkyView = 0 ; for i = 1:length(Hz) for j = 1:0.1:90 ; TotalSky = TotalSky + 1 ; if j > Hz(i) SkyView = SkyView + 1 ; end end end SVF = SkyView / TotalSky ; %Sun Beam (and weighted) View Factor calculation: TotalSunPath = 0 ; SunPathView = 0 ; TotalBeamEnergy = 0 ; BeamEnergyView = 0 ; for j=1:365 for i=1:1440 %Note: i represents the minute of the day and j represents the %day of the year. %The equation of time (E): B = ( j - 1 ) * 2 * pi / 365; E = 229.18 * ( 7.5e-5 + 1.868e-3 * cos(B) - 3.2077e-2 * sin(B)- 1.4615e-2 * cos(2*B) - 4.089e-2 * sin(2*B) ); %Convert to west longitude: if lon>0 longitude_standard=360-(time_zone+dst)*15; longitude_in_west_degree=360-lon; else longitude_standard= -(time_zone+dst)*15; longitude_in_west_degree= -lon; end %Difference in minutes between solar time and local time solartime_minus_localtime = 4 *( longitude_standard - longitude_in_west_degree )+ E; %Solar time (minute): solar_time= i + solartime_minus_localtime; if solar_time > 1440 solar_time = solar_time - 1440 ; end if solar_time <= 0 solar_time = 1440 - solar_time ; end %Decimal solar hour solar_hour = solar_time/60 ; %Solar angle calculation: %Hour angle (omega): omega = (solar_hour - 12) * 15 * pi / 180 ; %Declination angle (delta): delta = 23.45 * pi/180 * sin (2*pi*(j+284)/365); %Solar elevation angle (alpha_s): alpha_s = asin( cos(omega) * cos(delta) * cos(lat*pi/180) + sin(delta) * sin(lat*pi/180) ); alpha_s_degree = 180 * alpha_s/pi; % Solar azimuth angle (gamma_s): if lat >= 0 gamma_s = atan2( sin(omega) , cos(omega)*sin(lat*pi/180)-tan(delta)*cos(lat*pi/180) ) + pi ; else gamma_s = atan2( sin(omega) , cos(omega)*sin(lat*pi/180)-tan(delta)*cos(lat*pi/180) ) ; if gamma_s > 0 gamma_s = gamma_s - pi ; else gamma_s = gamma_s + pi ; end end gamma_s_degree = gamma_s * 180 / pi ; if alpha_s_degree > 90 theta_s = ( alpha_s_degree - 90 ) * pi / 180 ; else theta_s = ( 90 - alpha_s_degree ) * pi / 180 ; end if alpha_s > 0 RavR2 = 1.0001 + 0.034221 * cos(B) + 0.001280 * sin(B) + 0.000719 * cos(2*B) + 0.000077 * sin(2*B); %Note: RavR2 is the ratio beteween the mean sun-earth distance %(Rav) and the the actual sun-earth distance depending on the day of the year I0 = 1367 * RavR2 ; %Extraterrestrial radiation in W.m-2 AM = 1 / cos(theta_s) ; BNI = I0*0.7^(AM^0.678) ; TotalSunPath = TotalSunPath + 1 ; TotalBeamEnergy = TotalBeamEnergy + BNI ; ind = find( azi <= gamma_s ) ; if alpha_s_degree > Hz(ind(1)) SunPathView = SunPathView + 1 ; BeamEnergyView = BeamEnergyView + BNI ; end end end end SBVF = SunPathView / TotalSunPath ; SBWVF = BeamEnergyView / TotalBeamEnergy ; ViewF(n,m,1) = lat ; ViewF(n,m,2) = lon ; ViewF(n,m,3) = SVF ; ViewF(n,m,4) = SBVF ; ViewF(n,m,5) = SBWVF ; end end matlabpool('close'); toc
Merci par avance
Partager