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 :
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
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.
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 :

Error: The variable ViewF in a parfor cannot be classified.
Voici plus bas comment j'utilise ViewF:

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
Je ne comprends pas ce qui pose problème.

Merci par avance