Support for creating high altitude noctilucent clouds in Advanced Weather/ALS
[fg:toms-fgdata.git] / Nasal / local_weather / local_weather.nas
1
2 ########################################################
3 # routines to set up, transform and manage local weather
4 # Thorsten Renk, June 2011
5 # thermal model by Patrice Poly, April 2010
6 ########################################################
7
8 # function                      purpose
9 #
10 # calc_geo                      to compute the latitude to meter conversion
11 # calc_d_sq                     to compute a distance square in local Cartesian approximation
12 # effect_volume_loop            to check if the aircraft has entered an effect volume
13 # assemble_effect_array         to store the size of the effect volume array
14 # add_vectors                   to add two vectors in polar coordinates
15 # wind_altitude_interpolation   to interpolate aloft winds in altitude
16 # wind_interpolation            to interpolate aloft winds in altitude and position
17 # get_slowdown_fraction         to compute the effect of boundary layer wind slowdown
18 # interpolation_loop            to continuously interpolate weather parameters between stations 
19 # thermal_lift_start            to start the detailed thermal model
20 # thermal_lift_loop             to manage the detailed thermal lift model
21 # thermal_lift_stop             to end the detailed thermal lift model
22 # wave_lift_start               to start the detailed wave lift model
23 # wave_lift_loop                to manage the detailed wave lift model
24 # wave_lift_stop                to end the detailed wave lift model
25 # effect_volume_start           to manage parameters when an effect volume is entered
26 # effect_volume_stop            to manage parameters when an effect volume is left
27 # ts_factor                     (helper function for thermal lift model)
28 # tl_factor                     (helper function for thermal lift model)
29 # calcLift_max                  to calculate the maximal available thermal lift for given altitude
30 # calcLift                      to calculate the thermal lift at aircraft position
31 # calcWaveLift                  to calculate wave lift at aircraft position
32 # create_cloud_vec              to place a single cloud into an array to be written later
33 # clear_all                     to remove all clouds, effect volumes and weather stations and stop loops
34 # create_detailed_cumulus_cloud to place multiple cloudlets into a box based on a size parameter
35 # create_cumulonimbus_cloud     to place multiple cloudlets into a box 
36 # create_cumulonimbus_cloud_rain to place multiple cloudlets into a box and add a rain layer beneath
37 # create_cumosys                (wrapper to place a convective cloud system based on terrain coverage)
38 # cumulus_loop                  to place 25 Cumulus clouds each frame
39 # create_cumulus                to place a convective cloud system based on terrain coverage
40 # recreate_cumulus              to respawn convective clouds as part of the convective dynamics algorithm
41 # cumulus_exclusion_layer       to create a layer with 'holes' left for thunderstorm placement
42 # create_rise_clouds            to create a barrier cloud system
43 # create_streak                 to create a cloud streak
44 # create_hollow_layer           to create a cloud layer in a hollow cylinder (better for performance)
45 # create_cloudbox               to create a sophisticated cumulus cloud with different textures (experimental)
46 # terrain_presampling_start     to initialize terrain presampling
47 # terrain_presampling_loop      to sample 25 terrain points per frame
48 # terrain_presampling           to sample terrain elevation at a random point within specified area
49 # terrain_presampling_analysis  to analyze terrain presampling results
50 # wave_detection_loop           to detect if and where wave lift should be placed (currently unfinished)
51 # get_convective_altitude       to determine the altitude at which a Cumulus cloud is placed
52 # manage presampling            to take proper action when a presampling call has been finished
53 # set_wind_model_flag           to convert the wind model string into an integer flag
54 # set_texture_mix               to determine the texture mix between smooth and rough cloud appearance
55 # create_effect_volume          to create an effect volume
56 # set_weather_station           to specify a weather station for interpolation
57 # set_atmosphere_ipoint         to specify an interpolation point for visibility, haze and shading in the atmosphere
58 # set_wind_ipoint               to set an aloft wind interpolation point
59 # set_wind_ipoint_metar         to set a wind interpolation point from available ground METAR info where aloft is modelled
60 # showDialog                    to pop up a dialog window
61 # readFlags                     to read configuration flags from the property tree into Nasal variables at startup
62 # streak_wrapper                wrapper to execute streak from menu
63 # convection wrapper            wrapper to execute convective clouds from menu
64 # barrier wrapper               wrapper to execute barrier clouds from menu
65 # single_cloud_wrapper          wrapper to create single cloud from menu
66 # layer wrapper                 wrapper to create layer from menu
67 # box wrapper                   wrapper to create a cloudbox (experimental)
68 # set_aloft wrapper             wrapper to create aloft winds from menu
69 # set_tile                      to call a weather tile creation from menu
70 # startup                       to prepare the package at startup
71 # test                          to serve as a testbed for new functions
72
73 # object                        purpose
74
75 # weatherStation                to store info about weather conditions
76 # atmopshereIpoint              to store info about haze and light propagation in the atmosphere
77 # windIpoint                    to store an interpolation point of the windfield
78 # effectVolume                  to store effect volume info and provide methods to move and time-evolve effect volumes
79 # thermalLift                   to store thermal info and provide methods to move and time-evolve a thermal
80 # waveLift                      to store wave info 
81
82
83
84
85 ###################################
86 # geospatial helper functions
87 ###################################
88
89 var calc_geo = func(clat) {
90
91 lon_to_m  = math.cos(clat*math.pi/180.0) * lat_to_m;
92 m_to_lon = 1.0/lon_to_m;
93
94 weather_dynamics.lon_to_m = lon_to_m;
95 weather_dynamics.m_to_lon = m_to_lon;
96
97 }
98
99
100 var calc_d_sq = func (lat1, lon1, lat2, lon2) {
101
102 var x = (lat1 - lat2) * lat_to_m;
103 var y = (lon1 - lon2) * lon_to_m;
104
105 return (x*x + y*y);
106 }
107
108
109 ###################################
110 # effect volume management loop
111 ###################################
112
113 var effect_volume_loop = func (index, n_active) {
114
115
116 if (local_weather_running_flag == 0) {return;}
117
118 var n = 25;
119
120
121 var esize = n_effectVolumeArray;
122
123 var viewpos = geo.aircraft_position();
124 var active_counter = n_active;
125
126 var i_max = index + n;
127 if (i_max > esize) {i_max = esize;}
128
129 for (var i = index; i < i_max; i = i+1)
130         {
131         var e = effectVolumeArray[i];
132         
133         var flag = 0; # default assumption is that we're not in the volume
134         
135         var ealt_min = e.alt_low * ft_to_m;
136         var ealt_max = e.alt_high * ft_to_m;
137
138         
139         if ((viewpos.alt() > ealt_min) and (viewpos.alt() < ealt_max)) # we are in the correct alt range
140                 {
141                 # so we load geometry next
142                 
143                 var geometry = e.geometry;
144                 var elat = e.lat;
145                 var elon = e.lon;
146                 var rx = e.r1;
147
148                 if (geometry == 1) # we have a cylinder
149                         {
150                         var d_sq = calc_d_sq(viewpos.lat(), viewpos.lon(), elat, elon);
151                         if (d_sq < (rx*rx)) {flag =1;}
152                         }
153                 else if (geometry == 2) # we have an elliptic shape
154                         {
155                         # get orientation
156
157                         var ry = e.r2;
158                         var phi = e.phi;                
159
160                         phi = phi * math.pi/180.0;
161                                                 
162
163                         # first get unrotated coordinates 
164                         var xx = (viewpos.lon() - elon) * lon_to_m;
165                         var yy = (viewpos.lat() - elat) * lat_to_m;
166                         
167                         # then rotate to align with the shape
168                         var x = xx * math.cos(phi) - yy * math.sin(phi);
169                         var y = yy * math.cos(phi) + xx * math.sin(phi); 
170
171                         # then check elliptic condition
172                         if ((x*x)/(rx*rx) + (y*y)/(ry*ry) <1) {flag = 1;}
173                         }
174                 else if (geometry == 3) # we have a rectangular shape
175                         {
176                         # get orientation
177
178                         var ry = e.r2;
179                         var phi = e.phi;
180
181                         phi = phi * math.pi/180.0;
182                         # first get unrotated coordinates 
183                         var xx = (viewpos.lon() - elon) * lon_to_m;
184                         var yy = (viewpos.lat() - elat) * lat_to_m;
185                         # then rotate to align with the shape
186                         var x = xx * math.cos(phi) - yy * math.sin(phi);
187                         var y = yy * math.cos(phi) + xx * math.sin(phi); 
188                         # then check rectangle condition
189                         if ((x>-rx) and (x<rx) and (y>-ry) and (y<ry)) {flag = 1;}
190                         }
191                 } # end if altitude
192         
193         
194         # if flag ==1 at this point, we are inside the effect volume
195         # but we only need to take action on entering and leaving, so we check also active_flag
196         
197         # if (flag==1) {print("Inside volume");}
198         
199         var active_flag = e.active_flag;
200
201         if ((flag==1) and (active_flag ==0)) # we just entered the node
202                 {
203                 #print("Entered volume");               
204                 e.active_flag = 1;      
205                 effect_volume_start(e);
206                 }
207         else if ((flag==0) and (active_flag ==1)) # we left an active node
208                 {
209                 #print("Left volume!");
210                 e.active_flag = 0;
211                 effect_volume_stop(e);
212                 }
213         if (flag==1) {active_counter = active_counter + 1;} # we still count the active volumes
214         
215         } # end foreach
216
217 # at this point, all active effect counters should have been set to zero if we're outside all volumes
218 # however there seem to be rare configurations of overlapping volumes for which this doesn't happen
219 # therefore we zero them for redundancy here so that the interpolation loop can take over
220 # and set the properties correctly for outside
221
222
223 if (i == esize) # we check the number of actives and reset all counters
224         {
225         if (active_counter == 0)
226                 {
227                 var vNode = props.globals.getNode("local-weather/effect-volumes", 1);
228                 vNode.getChild("number-active-vis").setValue(0);
229                 vNode.getChild("number-active-snow").setValue(0);
230                 vNode.getChild("number-active-rain").setValue(0);
231                 vNode.getChild("number-active-lift").setValue(0);
232                 vNode.getChild("number-active-turb").setValue(0);
233                 vNode.getChild("number-active-sat").setValue(0);
234                 }
235         #print("n_active: ", active_counter);
236         active_counter = 0; i = 0;
237         }
238
239 # and we repeat the loop as long as the control flag is set
240
241
242 if (getprop(lw~"effect-loop-flag") ==1) {settimer( func {effect_volume_loop(i, active_counter); },0);}
243 }
244
245
246 ###################################
247 # assemble effect volume array
248 ###################################
249
250
251 var assemble_effect_array = func {
252
253 n_effectVolumeArray = size(effectVolumeArray);
254 }
255
256
257
258 ###################################
259 # vector addition
260 ###################################
261
262 var add_vectors = func (phi1, r1, phi2, r2) {
263
264 phi1 = phi1 * math.pi/180.0;
265 phi2 = phi2 * math.pi/180.0;
266
267 var x1 = r1 * math.sin(phi1);
268 var x2 = r2 * math.sin(phi2);
269
270 var y1 = r1 * math.cos(phi1);
271 var y2 = r2 * math.cos(phi2);
272
273 var x = x1+x2;
274 var y = y1+y2;
275
276 var phi = math.atan2(x,y) * 180.0/math.pi;
277 var r = math.sqrt(x*x + y*y);
278
279 var vec = [];
280
281 append(vec, phi);
282 append(vec,r);
283
284 return vec;
285 }
286
287
288 ###################################
289 # windfield altitude interpolation
290 ###################################
291
292
293 var wind_altitude_interpolation = func (altitude, w) {
294
295 if (altitude < wind_altitude_array[0]) {var alt_wind = wind_altitude_array[0];}
296 else if (altitude > wind_altitude_array[8]) {var alt_wind = 0.99* wind_altitude_array[8];}
297 else {var alt_wind = altitude;}
298
299 for (var i = 0; i<9; i=i+1)
300         {if (alt_wind < wind_altitude_array[i]) {break;}}
301         
302
303 #var altNodeMin = w.getChild("altitude",i-1);
304 #var altNodeMax = w.getChild("altitude",i);     
305
306 #var vmin = altNodeMin.getNode("windspeed-kt").getValue();
307 #var vmax = altNodeMax.getNode("windspeed-kt").getValue();
308
309 var vmin = w.alt[i-1].v;
310 var vmax = w.alt[i].v;
311
312 #var dir_min = altNodeMin.getNode("wind-from-heading-deg").getValue();
313 #var dir_max = altNodeMax.getNode("wind-from-heading-deg").getValue();
314
315 var dir_min = w.alt[i-1].d;
316 var dir_max = w.alt[i].d;
317
318 var f = (alt_wind - wind_altitude_array[i-1])/(wind_altitude_array[i] - wind_altitude_array[i-1]);
319
320 var res = add_vectors(dir_min, (1-f) * vmin, dir_max, f * vmax);
321
322 return res;
323 }
324
325
326 ###################################
327 # windfield spatial interpolation
328 ###################################
329
330 var wind_interpolation = func (lat, lon, alt) {
331
332 var sum_norm = 0;
333 var sum_wind = [0,0];
334
335 var wsize = size(windIpointArray);
336         
337 for (var i = 0; i < wsize; i=i+1) {
338         
339         
340         var w = windIpointArray[i];
341
342         var wpos = geo.Coord.new();
343         wpos.set_latlon(w.lat,w.lon,1000.0);
344
345         var ppos = geo.Coord.new();
346         ppos.set_latlon(lat,lon,1000.0);
347
348         var d = ppos.distance_to(wpos);
349         if (d <100.0) {d = 100.0;} # to prevent singularity at zero
350
351         sum_norm = sum_norm + (1./d) * w.weight;
352
353         var res = wind_altitude_interpolation(alt,w);
354         
355         sum_wind = add_vectors(sum_wind[0], sum_wind[1], res[0], (res[1]/d) * w.weight);        
356
357         # gradually fade in the interpolation weight of newly added points to
358         # avoid sudden jumps
359
360         if (w.weight < 1.0) {w.weight = w.weight + 0.02;}
361
362         }
363
364 sum_wind[1] = sum_wind[1] /sum_norm;
365
366 return sum_wind;
367 }
368
369
370 ###################################
371 # boundary layer computations
372 ###################################
373
374
375 var get_slowdown_fraction = func {
376
377 var tile_index = getprop(lw~"tiles/tile[4]/tile-index");
378 var altitude_agl = getprop("/position/altitude-agl-ft");
379 var altitude = getprop("/position/altitude-ft");
380
381
382
383 if (presampling_flag == 0)
384         {
385         var base_layer_thickness = 600.0;       
386         var f_slow = 1.0/3.0;
387         }
388 else 
389         {
390         var alt_median = alt_50_array[tile_index - 1];
391         var alt_difference = alt_median - (altitude - altitude_agl);
392         var base_layer_thickness = 150.0;       
393
394         # get the boundary layer size dependent on terrain altitude above terrain median
395
396         if (alt_difference > 0.0) # we're low and the boundary layer grows
397                 {var boundary_alt = base_layer_thickness + 0.3 * alt_difference;}
398         else # the boundary layer shrinks
399                 {var boundary_alt = base_layer_thickness + 0.1 * alt_difference;}
400
401         if (boundary_alt < 50.0){boundary_alt = 50.0;}
402         if (boundary_alt > 3000.0) {boundary_alt = 3000.0;}
403
404         # get the boundary effect as a function of bounday layer size
405         
406         var f_slow = 1.0 - (0.2 + 0.17 * math.ln(boundary_alt/base_layer_thickness));
407         }
408
409 if (debug_output_flag == 1)
410         {
411         #print("Boundary layer thickness: ",base_layer_thickness);
412         #print("Boundary layer slowdown: ", f_slow);
413         }
414 return f_slow;
415 }
416
417
418 ###################################
419 # interpolation management loop
420 ###################################
421
422 var interpolation_loop = func {
423
424 if (local_weather_running_flag == 0) {return;}
425
426 var viewpos = geo.aircraft_position();
427
428
429
430 #var vis_before = getprop(lwi~"visibility-m");
431 var vis_before = interpolated_conditions.visibility_m;
432
433 # if applicable, do some work for fps sampling
434
435 if (fps_control_flag == 1)
436         {
437         fps_samples = fps_samples +1;
438         fps_sum = fps_sum + getprop("/sim/frame-rate");
439         }
440
441
442 # determine at which distance we no longer keep an interpolation point, needs to be larger for METAR since points are more scarce
443
444 if (metar_flag == 1)
445         {var distance_to_unload = 250000.0;}
446 else    
447         {var distance_to_unload = 120000.0;}    
448
449 # if we can set environment without a reset, the loop can run a bit faster for smoother interpolation
450 # so determine the suitable timing
451
452
453 var interpolation_loop_time = 0.1; 
454 var vlimit = 1.01;
455
456
457 # get an inverse distance weighted average from all defined weather stations
458
459 var sum_alt = 0.0;
460 var sum_vis = 0.0;
461 var sum_T = 0.0;
462 var sum_p = 0.0;
463 var sum_D = 0.0;
464 var sum_norm = 0.0;
465
466 var n_stations = size(weatherStationArray);
467
468 for (var i = 0; i < n_stations; i=i+1) {
469         
470         var s = weatherStationArray[i];
471         
472
473         var stpos = geo.Coord.new();
474         stpos.set_latlon(s.lat,s.lon,0.0);
475
476         var d = viewpos.distance_to(stpos);
477         if (d <100.0) {d = 100.0;} # to prevent singularity at zero
478
479         sum_norm = sum_norm + 1./d * s.weight;
480         
481         sum_alt = sum_alt + (s.alt/d) * s.weight;
482         sum_vis = sum_vis + (s.vis/d) * s.weight;
483         sum_T = sum_T + (s.T/d) * s.weight;
484         sum_D = sum_D + (s.D/d) * s.weight;
485         sum_p = sum_p + (s.p/d) * s.weight;
486
487         # gradually fade in the interpolation weight of newly added stations to
488         # avoid sudden jumps
489
490         if (s.weight < 1.0) {s.weight = s.weight + 0.02;}
491
492         # automatically delete stations out of range
493         # take care not to unload if weird values appear for a moment
494         # never unload if only one station left
495         if ((d > distance_to_unload) and (d < (distance_to_unload + 20000.0)) and (n_stations > 1)) 
496                 {
497                 if (debug_output_flag == 1) 
498                         {print("Distance to weather station ", d, " m, unloading ...", i);}
499                 weatherStationArray = weather_tile_management.delete_from_vector(weatherStationArray,i);
500                 i = i-1; n_stations = n_stations -1;
501                 }
502         }
503
504 setprop(lwi~"station-number", i+1);
505
506
507 var ialt = sum_alt/sum_norm;
508 var vis = sum_vis/sum_norm;
509 var p = sum_p/sum_norm;
510 var D = sum_D/sum_norm + temperature_offset;
511 var T = sum_T/sum_norm + temperature_offset;
512
513
514
515 # get an inverse distance weighted average from all defined atmospheric condition points
516
517 sum_norm = 0.0;
518 var sum_vis_aloft = 0.0;
519 var sum_vis_alt1 = 0.0;
520 var sum_vis_ovcst = 0.0;
521 var sum_ovcst = 0.0;
522 var sum_ovcst_alt_low = 0.0;
523 var sum_ovcst_alt_high = 0.0;
524 var sum_scatt = 0.0;
525 var sum_scatt_alt_low = 0.0;
526 var sum_scatt_alt_high = 0.0;
527
528 var n_iPoints = size(atmosphereIpointArray);
529
530 for (var i = 0; i < n_iPoints; i=i+1) {
531         
532         var a = atmosphereIpointArray[i];
533         
534
535         var apos = geo.Coord.new();
536         apos.set_latlon(a.lat,a.lon,0.0);
537
538         var d = viewpos.distance_to(apos);
539         if (d <100.0) {d = 100.0;} # to prevent singularity at zero
540
541         sum_norm = sum_norm + 1./d * a.weight;
542         sum_vis_aloft = sum_vis_aloft + (a.vis_aloft/d) * a.weight;
543         sum_vis_alt1 = sum_vis_alt1 + (a.vis_alt1/d) * a.weight;
544         sum_vis_ovcst = sum_vis_ovcst + (a.vis_ovcst/d) * a.weight;     
545         sum_ovcst = sum_ovcst + (a.ovcst/d) * a.weight;
546         sum_ovcst_alt_low = sum_ovcst_alt_low + (a.ovcst_alt_low/d) * a.weight;
547         sum_ovcst_alt_high = sum_ovcst_alt_high + (a.ovcst_alt_high/d) * a.weight;
548         sum_scatt = sum_scatt + (a.scatt/d) * a.weight;
549         sum_scatt_alt_low = sum_scatt_alt_low + (a.scatt_alt_low/d) * a.weight;
550         sum_scatt_alt_high = sum_scatt_alt_high + (a.scatt_alt_high/d) * a.weight;
551
552         # gradually fade in the interpolation weight of newly added stations to
553         # avoid sudden jumps
554
555         if (a.weight < 1.0) {a.weight = a.weight + 0.02;}
556
557         # automatically delete stations out of range
558         # take care not to unload if weird values appear for a moment
559         # never unload if only one station left
560         if ((d > distance_to_unload) and (d < (distance_to_unload + 20000.0)) and (n_iPoints > 1)) 
561                 {
562                 if (debug_output_flag == 1) 
563                         {print("Distance to atmosphere interpolation point ", d, " m, unloading ...", i);}
564                 atmosphereIpointArray = weather_tile_management.delete_from_vector(atmosphereIpointArray,i);
565                 i = i-1; n_iPoints = n_iPoints -1;
566                 }
567         }
568
569 setprop(lwi~"atmosphere-ipoint-number", i+1);
570
571
572
573 var vis_aloft = sum_vis_aloft/sum_norm;
574 var vis_alt1 = sum_vis_alt1/sum_norm;
575 var vis_ovcst = sum_vis_ovcst/sum_norm;
576 var ovcst_max = sum_ovcst/sum_norm;
577 var ovcst_alt_low = sum_ovcst_alt_low/sum_norm;
578 var ovcst_alt_high = sum_ovcst_alt_high/sum_norm;
579 var scatt_max = sum_scatt/sum_norm;
580 var scatt_alt_low = sum_scatt_alt_low/sum_norm;
581 var scatt_alt_high = sum_scatt_alt_high/sum_norm;
582
583
584
585
586 # altitude model for visibility - increase above the lowest inversion layer to simulate ground haze
587
588 vis = vis * ground_haze_factor;
589
590 var altitude = getprop("position/altitude-ft");
591 # var current_mean_terrain_elevation = ialt;
592
593 var alt1 = vis_alt1;
594 var alt2 = alt1 + 1500.0;
595
596
597 setprop("/environment/ground-visibility-m",vis);
598 setprop("/environment/ground-haze-thickness-m",alt2 * ft_to_m);
599
600 # compute the visibility gradients
601
602 if (realistic_visibility_flag == 1)
603         {
604         vis_aloft = vis_aloft * 2.0;
605         vis_ovcst = vis_ovcst * 3.0;
606         }
607
608 var inc1 = 0.0 * (vis_aloft - vis)/(vis_alt1 - ialt);
609 var inc2 = 0.9 * (vis_aloft - vis)/1500.0;
610 var inc3 = (vis_ovcst - vis_aloft)/(ovcst_alt_high - vis_alt1+1500);
611 var inc4 = 0.5;
612
613
614 if (realistic_visibility_flag == 1)
615         {inc4 = inc4 * 3.0;}
616
617 # compute the visibility
618
619 if (altitude < alt1)
620         {vis = vis + inc1 * altitude;}
621 else if (altitude < alt2)
622         {
623         vis = vis + inc1 * alt1 + inc2 * (altitude - alt1); 
624         }
625 else if (altitude < ovcst_alt_high)
626         {
627         vis = vis + inc1 * alt1 + inc2 * (alt2-alt1)  + inc3 * (altitude - alt2);
628         }
629 else if (altitude > ovcst_alt_high)
630         {
631         vis = vis + inc1 * alt1 + inc2 * (alt2-alt1)  + inc3 * (ovcst_alt_high - alt2) + inc4 * (altitude - ovcst_alt_high);
632         }
633
634 # limit visibility (otherwise memory consumption may be very bad...)
635
636 if (vis > max_vis_range)
637         {vis = max_vis_range;}
638
639         
640 # determine scattering shader parameters if scattering shader is on
641
642 if (scattering_shader_flag == 1) 
643         {
644         
645         # values to be used with new exposure filter
646         var rayleigh = 0.0003;
647         var mie = 0.005;
648         var density = 0.3;
649
650         var vis_factor = (vis - 30000.0)/90000.0;
651         if (vis_factor < 0.0) {vis_factor = 0.0;}
652         if (vis_factor > 1.0) {vis_factor = 1.0;}
653  
654
655         if (altitude < 36000.0) 
656                 {
657                 rayleigh = 0.0003 - 0.0001 * vis_factor;
658                 mie = 0.005 - vis_factor * 0.002; 
659                 }
660         else if (altitude < 85000.0)
661                 {
662                 rayleigh = (0.0003 - 0.0001 * vis_factor)  - (altitude-36000.0)/49000.0 * 0.0001;
663                 mie = 0.005 - vis_factor * 0.002 - (altitude-36000.0)/49000.0 * 0.002;
664                 }
665         else 
666                 {rayleigh = 0.0002 - 0.0001 * vis_factor; mie = 0.003 - vis_factor * 0.002;}
667
668        # now the pollution factor
669    
670         if (altitude < alt1)
671                 {
672                 rayleigh = rayleigh +0.0003 * air_pollution_norm + 0.0004 * air_pollution_norm * (1.0 - (altitude/alt1) * (altitude/alt1));
673                 density = density + 0.05 * air_pollution_norm + 0.05 * air_pollution_norm * (1.0 - (altitude/alt1) * (altitude/alt1));
674                 }
675         else
676                 {
677                 rayleigh = rayleigh + 0.0003 * air_pollution_norm;
678                 density = density + 0.05 * air_pollution_norm;
679                 }
680
681
682         }
683
684
685 # compute the horizon shading
686
687 if (altitude < scatt_alt_low)
688         {
689         var scatt = scatt_max;
690         }
691 else if (altitude < scatt_alt_high)
692         {
693         var scatt = scatt_max + (0.95 - scatt_max) * (altitude - scatt_alt_low)/(scatt_alt_high-scatt_alt_low);
694         }
695 else
696         {var scatt = 0.95;}
697
698
699 # compute  the cloud layer self shading correction
700
701 var sun_angle = 1.57079632675 - getprop("/sim/time/sun-angle-rad");
702 var cloud_layer_shading = 1.0 - (0.8*(1.0 - scatt_max) *  math.pow(math.cos(sun_angle),100.0));
703
704 # compute the overcast haze
705
706 if (altitude < ovcst_alt_low)
707         {
708         var ovcst = ovcst_max;
709         }       
710 else if (altitude < ovcst_alt_high)
711         {
712         var ovcst = ovcst_max - ovcst_max * (altitude - ovcst_alt_low)/(ovcst_alt_high-ovcst_alt_low);
713         }
714 else
715         {var ovcst = 0.0;}
716
717
718 # compute base turbulence
719
720 var base_turbulence = 0.0;
721
722 if (altitude < alt1)
723         {
724         base_turbulence = lowest_layer_turbulence;
725         }
726
727
728
729 # limit relative changes of the visibility, will make for gradual transitions
730
731 if (vis/vis_before > vlimit)
732         {vis = vlimit * vis_before;}
733 else if (vis/vis_before < (2.0-vlimit))
734         {vis = (2.0-vlimit) * vis_before;}
735
736
737
738
739 # write all properties into the weather interpolation record 
740
741 setprop(lwi~"mean-terrain-altitude-ft",ialt);
742
743
744 if (vis > 0.0) interpolated_conditions.visibility_m = vis;
745 interpolated_conditions.temperature_degc = T;
746 interpolated_conditions.dewpoint_degc = D;
747 if (p>10.0) interpolated_conditions.pressure_sea_level_inhg = p;
748
749
750
751 if (scattering_shader_flag == 1)
752         {
753         local_weather.setSkydomeShader(rayleigh, mie, density);
754         setprop("/environment/cloud-self-shading", cloud_layer_shading);
755         }
756
757 local_weather.setScattering(scatt);
758 local_weather.setOvercast(ovcst);
759
760
761         
762
763 # now check if an effect volume writes the property and set only if not
764 # but set visibility if interpolated is smaller than effect-specified
765
766 var flag = getprop("local-weather/effect-volumes/number-active-vis");
767
768 if ((flag ==0) and (vis > 0.0) and (getprop(lw~"lift-loop-flag") == 0) and (compat_layer.smooth_visibility_loop_flag == 0))
769         {
770         compat_layer.setVisibility(vis);
771         }
772 else if ((getprop("/local-weather/current/visibility-m") > vis) and (compat_layer.smooth_visibility_loop_flag == 0))
773         {
774         compat_layer.setVisibility(vis);
775         }
776
777
778
779
780
781 flag = getprop("local-weather/effect-volumes/number-active-lift");
782
783 if (flag ==0) 
784         {
785         #setprop(lw~"current/thermal-lift",0.0);
786         }
787
788 # no need to check for these, as they are not modelled in effect volumes
789
790 compat_layer.setTemperature(T);
791 compat_layer.setDewpoint(D);
792 if (p>0.0) {compat_layer.setPressure(p);}
793
794
795 # now determine the local wind 
796
797
798 var tile_index = getprop(lw~"tiles/tile[4]/tile-index");
799
800 if (wind_model_flag ==1) # constant
801         {
802         var winddir = weather_dynamics.tile_wind_direction[0];
803         var windspeed = weather_dynamics.tile_wind_speed[0];
804
805         wind.cloudlayer = [winddir,windspeed];
806
807         }
808 else if (wind_model_flag ==2) # constant in tile
809         {
810         var winddir = weather_dynamics.tile_wind_direction[tile_index-1];
811         var windspeed = weather_dynamics.tile_wind_speed[tile_index-1];
812
813         wind.cloudlayer = [winddir, windspeed];
814
815         }       
816 else if (wind_model_flag ==3) # aloft interpolated, constant in tiles
817         {
818         var w = windIpointArray[0];
819         var res = wind_altitude_interpolation(altitude,w);
820         var winddir = res[0];
821         var windspeed = res[1];
822
823         wind.cloudlayer = wind_altitude_interpolation(0.0,w);
824
825         }
826 else if (wind_model_flag == 5) # aloft waypoint interpolated
827         {
828         var res = wind_interpolation(viewpos.lat(), viewpos.lon(), altitude);   
829
830         var winddir = res[0];
831         var windspeed = res[1];
832
833         wind.cloudlayer = wind_interpolation(viewpos.lat(), viewpos.lon(), 0.0);        
834         }
835
836
837 wind.surface = [wind.cloudlayer[0], wind.cloudlayer[1] * get_slowdown_fraction()];
838
839 # now do the boundary layer computations
840
841 var altitude_agl = getprop("/position/altitude-agl-ft");
842
843 if (altitude_agl < 50.0)
844         {
845         base_turbulence = base_turbulence * altitude_agl/50.0;
846         }
847
848
849 if (presampling_flag == 0)
850         {
851         var boundary_alt = 600.0;
852         var windspeed_ground = windspeed/3.0;
853         
854         var f_min = 2.0/3.0;
855
856         if (altitude_agl < boundary_alt)
857                 {var windspeed_current = windspeed_ground + 2.0 * windspeed_ground * (altitude_agl/boundary_alt);}
858         else 
859                 {var windspeed_current = windspeed;}
860         }
861 else 
862         {
863         var alt_median = alt_50_array[tile_index - 1];
864         var alt_difference = alt_median - (altitude - altitude_agl);
865         var base_layer_thickness = 150.0;       
866
867         # get the boundary layer size dependent on terrain altitude above terrain median
868
869         if (alt_difference > 0.0) # we're low and the boundary layer grows
870                 {var boundary_alt = base_layer_thickness + 0.3 * alt_difference;}
871         else # the boundary layer shrinks
872                 {var boundary_alt = base_layer_thickness + 0.1 * alt_difference;}
873
874         if (boundary_alt < 50.0){boundary_alt = 50.0;}
875         if (boundary_alt > 3000.0) {boundary_alt = 3000.0;}
876
877         # get the boundary effect as a function of bounday layer size
878         
879         var f_min = 0.2 + 0.17 * math.ln(boundary_alt/base_layer_thickness);
880
881
882         if (altitude_agl < boundary_alt)
883                 {
884                 var windspeed_current = (1-f_min) * windspeed + f_min * windspeed * (altitude_agl/boundary_alt);
885                 }
886         else 
887                 {var windspeed_current = windspeed;}
888
889         }
890
891
892 var windspeed_ground = (1.0-f_min) * windspeed;
893
894
895 # set the wind hash before gusts, it represents mean wind
896
897 wind.current = [winddir,windspeed_current];
898
899
900
901 # determine gusts and turbulence in the bounday layer
902
903 var gust_frequency = getprop(lw~"tmp/gust-frequency-hz");
904
905
906
907
908 if (gust_frequency > 0.0)
909         {
910         var gust_relative_strength = getprop(lw~"tmp/gust-relative-strength");
911         var gust_angvar = getprop(lw~"tmp/gust-angular-variation-deg");
912         
913         # var winddir_last = getprop(lwi~"wind-from-heading-deg");
914         var winddir_last = interpolated_conditions.wind_from_heading_deg;
915         
916         var alt_scaling_factor = 1.2 * windspeed / 10.0;
917         if (alt_scaling_factor < 1.0) {alt_scaling_factor = 1.0;}
918
919         # expected mean number of gusts in time interval (should be < 1)
920         var p_gust = gust_frequency * interpolation_loop_time;
921
922         winddir_change = 0.0;
923
924         if (rand() < p_gust) # we change the offsets for windspeed and direction
925                 {
926                 var alt_fact = 1.0 - altitude_agl/(boundary_alt * alt_scaling_factor);
927                 if (alt_fact < 0.0) {alt_fact = 0.0};
928                 windspeed_multiplier =  (1.0 + ((rand()) * gust_relative_strength * alt_fact));
929                 winddir_change = alt_fact * (1.0 - 2.0 * rand()) * gust_angvar;
930                 winddir_change = winddir_change * 0.2; # Markov chain parameter, max. change per frame is 1/5 
931                 
932                 # if the Markov chain reaches the boundary, reflect
933
934                 #print("Winddir: ", winddir, " winddir_last: ", winddir_last, " winddir_change: ", winddir_change);
935                 if (weather_tile_management.relangle(winddir_last + winddir_change, winddir) > gust_angvar)
936                         {winddir_change = -winddir_change;}
937                 
938                 }
939         windspeed_current = windspeed_current *  windspeed_multiplier;
940         winddir = winddir_last + winddir_change;
941         }
942
943         
944
945
946
947 compat_layer.setWindSmoothly(winddir, windspeed_current);
948
949 interpolated_conditions.wind_from_heading_deg = winddir;
950 interpolated_conditions.windspeed_kt = windspeed_current;
951
952 # hack to get access to the water shader
953
954 setprop("/environment/config/boundary/entry[0]/wind-from-heading-deg",winddir);
955 setprop("/environment/config/boundary/entry[0]/wind-speed-kt",windspeed_ground);
956
957 # end hack
958
959
960
961
962 # set turbulence
963 flag = getprop("local-weather/effect-volumes/number-active-turb");
964
965 var wind_enhancement_factor = windspeed_current/15.0;
966 if (wind_enhancement_factor > 1.5) {wind_enhancement_factor = 1.5;}
967
968 if ((flag ==0))
969         {compat_layer.setTurbulence(base_turbulence * wind_enhancement_factor);}
970
971 # set scattering on the ground - this doesn't affect fog but is diffuse and specular light reduction
972 # so it is stronger than normal scattering
973
974 var scatt_ground = (scatt_max - 0.4)/0.6;
975 if (scatt_ground < 0.0) {scatt_ground = 0.0;}
976
977 setprop("/environment/surface/scattering", scatt_ground);
978
979 if (getprop(lw~"interpolation-loop-flag") ==1) {settimer(interpolation_loop, 0.0);}
980
981 }
982
983 ###################################
984 # thermal lift loop startup
985 ###################################
986
987 var thermal_lift_start = func (ev) {
988
989
990 # if another lift loop is already running, do nothing
991 if (getprop(lw~"lift-loop-flag") == 1) {return;} 
992
993 # copy the properties from effect volume to the lift object
994
995 l = thermalLift.new(ev.lat, ev.lon, ev.radius, ev.height, ev.cn, ev.sh, ev.max_lift, ev.f_lift_radius);
996
997 l.index = ev.index;
998
999 if (dynamics_flag == 1)
1000         {
1001         l.timestamp = weather_dynamics.time_lw;
1002         if (dynamical_convection_flag == 1)
1003                 {
1004                 l.flt = ev.flt;
1005                 l.evolution_timestamp = ev.evolution_timestamp;
1006                 }
1007         }
1008
1009
1010
1011 thermal = l;
1012
1013 if (debug_output_flag == 1)
1014         {
1015         print("Entering thermal lift...");
1016         print("strength: ", thermal.max_lift, " radius: ", thermal.radius);
1017         if (dynamical_convection_flag ==1)
1018                 {print("fractional lifetime: ", thermal.flt);}
1019
1020         }
1021
1022 # and start the lift loop, unless another one is already running
1023 # so we block overlapping calls
1024
1025
1026 setprop(lw~"lift-loop-flag",1); 
1027 settimer(thermal_lift_loop,0);
1028
1029 }
1030
1031 ###################################
1032 # thermal lift loop
1033 ###################################
1034
1035 var thermal_lift_loop = func {
1036
1037 if (local_weather_running_flag == 0) {return;}
1038
1039 var apos = geo.aircraft_position();
1040
1041 var tlat = thermal.lat;
1042 var tlon = thermal.lon;
1043
1044 var tpos = geo.Coord.new();
1045 tpos.set_latlon(tlat,tlon,0.0);
1046
1047 var d = apos.distance_to(tpos);
1048 var alt = getprop("position/altitude-ft");
1049
1050 if (dynamical_convection_flag == 1)
1051         {var flt = thermal.flt;}
1052 else
1053         {var flt = 0.5;}
1054
1055 var lift = calcLift(d, alt, thermal.radius, thermal.height, thermal.cn, thermal.sh, thermal.max_lift, thermal.f_lift_radius, flt);
1056
1057 if (getprop(lw~"wave-loop-flag") ==1) 
1058         {
1059         lift = lift + getprop(lw~"current/wave-lift");
1060         }
1061
1062 # compute a reduction in visibility when entering the cloudbase
1063
1064 #var vis = getprop(lw~"interpolation/visibility-m");
1065
1066 var vis = interpolated_conditions.visibility_m;
1067
1068 if (alt > 0.9 * thermal.height)
1069         {
1070         var visibility_reduction = math.pow((alt - 0.9 * thermal.height)/(0.2 * thermal.height),0.1);
1071         visibility_reduction = visibility_reduction * (1.0 - math.pow(d/(0.8*thermal.radius),14));
1072
1073         if (visibility_reduction > 1.0) {visibility_reduction = 1.0;} # this shouldn't ever happen
1074         if (visibility_reduction < 0.0) {visibility_reduction = 0.0;} 
1075         vis = vis * (1.0 - 0.98 * visibility_reduction);
1076
1077         }
1078
1079 setprop(lw~"current/visibility-m",vis);
1080 compat_layer.setVisibility(vis);
1081
1082
1083
1084
1085 setprop(lw~"current/thermal-lift",lift);
1086 compat_layer.setLift(lift);
1087
1088 # if dynamics is on, move the thermal and occasionally compute altitude and age
1089
1090 if (dynamics_flag == 1)
1091         {
1092         thermal.move();
1093         
1094         if ((rand() < 0.01) and (presampling_flag == 1)) # check every 100 frames
1095                 {
1096                 if (dynamical_convection_flag == 1) 
1097                         {
1098                         thermal.correct_altitude_and_age();
1099                         if (thermal.flt > 1.1)
1100                                 {thermal_lift_stop();}
1101                         }
1102                 else    
1103                         {       
1104                         thermal.correct_altitude();
1105                         }
1106                 }       
1107         }
1108
1109
1110 if (getprop(lw~"lift-loop-flag") ==1) {settimer(thermal_lift_loop, 0);}
1111 }
1112
1113
1114
1115
1116
1117 ###################################
1118 # thermal lift loop stop
1119 ###################################
1120
1121 var thermal_lift_stop = func {
1122
1123 setprop(lw~"lift-loop-flag",0);
1124 setprop(lw~"current/thermal-lift",0.0);
1125 compat_layer.setLift(0.0);
1126
1127 if (debug_output_flag == 1)
1128         {
1129         print("Leaving thermal lift...");
1130         }
1131
1132 }
1133
1134
1135 ###################################
1136 # wave lift loop startup
1137 ###################################
1138
1139 var wave_lift_start = func (ev) {
1140
1141 # copy the properties from effect volume to the wave object
1142
1143
1144 w = waveLift.new (ev.lat, ev.lon, ev.r1, ev.r2, ev.phi, ev.height, ev.max_lift);
1145 w.index = ev.index;
1146 wave = w;
1147
1148 # and start the lift loop, unless another one is already running
1149 # so we block overlapping calls
1150
1151 if (getprop(lw~"wave-loop-flag") == 0) 
1152 {setprop(lw~"wave-loop-flag",1); settimer(wave_lift_loop,0);}
1153
1154 }
1155
1156 ###################################
1157 # wave lift loop
1158 ###################################
1159
1160 var wave_lift_loop = func {
1161
1162 if (local_weather_running_flag == 0) {return;}
1163
1164 var lat = getprop("position/latitude-deg");
1165 var lon = getprop("position/longitude-deg");
1166 var alt = getprop("position/altitude-ft");
1167
1168
1169 var phi = wave.phi * math.pi/180.0;
1170
1171 var xx = (lon - wave.lon) * lon_to_m;
1172 var yy = (lat - wave.lat) * lat_to_m;
1173
1174 var x = xx * math.cos(phi) - yy * math.sin(phi);
1175 var y = yy * math.cos(phi) + xx * math.sin(phi); 
1176
1177 var lift = calcWaveLift(x,y,alt);
1178
1179 # check if we are in a thermal, if so set wave lift and let the thermal lift loop add that
1180
1181 if (getprop(lw~"lift-loop-flag") ==1)
1182         {
1183         setprop(lw~"current/wave-lift",lift);
1184         }
1185 else
1186         {
1187         setprop(lw~"current/thermal-lift",lift);
1188         }
1189
1190 if (getprop(lw~"wave-loop-flag") ==1) {settimer(wave_lift_loop, 0);}
1191 }
1192
1193
1194
1195
1196 ###################################
1197 # wave lift loop stop
1198 ###################################
1199
1200 var wave_lift_stop = func {
1201
1202 setprop(lw~"wave-loop-flag",0);
1203 setprop(lw~"current/thermal-lift",0.0);
1204 }
1205
1206
1207
1208 ####################################
1209 # action taken when in effect volume
1210 ####################################
1211
1212 var effect_volume_start = func (ev) {
1213
1214 var cNode = props.globals.getNode(lw~"current");
1215
1216
1217 if (ev.vis_flag ==1)
1218         {
1219         # first store the current setting in case we need to restore on leaving 
1220         
1221         var vis = ev.vis;
1222         ev.vis_r = cNode.getNode("visibility-m").getValue();
1223
1224         # then set the new value in current and execute change
1225         cNode.getNode("visibility-m").setValue(vis);
1226         #compat_layer.setVisibility(vis);
1227         print(vis);
1228         compat_layer.setVisibilitySmoothly(vis);
1229
1230         # then count the number of active volumes on entry (we need that to determine
1231         # what to do on exit)
1232         ev.n_entry_vis = getprop(lw~"effect-volumes/number-active-vis");
1233
1234         # and add to the counter
1235         setprop(lw~"effect-volumes/number-active-vis",getprop(lw~"effect-volumes/number-active-vis")+1);
1236         }
1237
1238 if (ev.rain_flag == 1)
1239         {
1240         var rain = ev.rain;
1241         #print("Setting rain to:", rain);
1242         ev.rain_r = cNode.getNode("rain-norm").getValue();
1243         cNode.getNode("rain-norm").setValue(rain);
1244         compat_layer.setRain(rain);
1245         ev.n_entry_rain = getprop(lw~"effect-volumes/number-active-rain");
1246         setprop(lw~"effect-volumes/number-active-rain",getprop(lw~"effect-volumes/number-active-rain")+1);
1247         }
1248 if (ev.snow_flag == 1)
1249         {
1250         var snow = ev.snow;
1251         ev.snow_r = cNode.getNode("snow-norm").getValue();
1252         cNode.getNode("snow-norm").setValue(snow);
1253         compat_layer.setSnow(snow);
1254         ev.n_entry_snow = getprop(lw~"effect-volumes/number-active-snow");
1255         setprop(lw~"effect-volumes/number-active-snow",getprop(lw~"effect-volumes/number-active-snow")+1);
1256         }
1257 if (ev.turb_flag == 1)
1258         {
1259         var turbulence = ev.turb;
1260         ev.turb_r = cNode.getNode("turbulence").getValue();
1261         cNode.getNode("turbulence").setValue(turbulence);
1262         compat_layer.setTurbulence(turbulence);
1263         ev.n_entry_turb = getprop(lw~"effect-volumes/number-active-turb");
1264         setprop(lw~"effect-volumes/number-active-turb",getprop(lw~"effect-volumes/number-active-turb")+1);
1265         }
1266 if (ev.sat_flag == 1)
1267         {
1268         var saturation = ev.sat;
1269         ev.sat_r = getprop("/rendering/scene/saturation");
1270         compat_layer.setLightSmoothly(saturation);
1271         ev.n_entry_sat = getprop(lw~"effect-volumes/number-active-sat");
1272         setprop(lw~"effect-volumes/number-active-sat",getprop(lw~"effect-volumes/number-active-sat")+1);
1273         }
1274
1275 if (ev.lift_flag == 1)
1276         {
1277         var lift = ev.lift;
1278         ev.lift_r = cNode.getNode("thermal-lift").getValue();
1279         cNode.getNode("thermal-lift").setValue(lift);
1280         compat_layer.setLift(lift);
1281         ev.n_entry_lift = getprop(lw~"effect-volumes/number-active-lift");      
1282         setprop(lw~"effect-volumes/number-active-lift",getprop(lw~"effect-volumes/number-active-lift")+1);
1283         }
1284 else if (ev.lift_flag == 2)
1285         {
1286         ev.lift_r = cNode.getNode("thermal-lift").getValue();
1287         ev.n_entry_lift = getprop(lw~"effect-volumes/number-active-lift");      
1288         setprop(lw~"effect-volumes/number-active-lift",getprop(lw~"effect-volumes/number-active-lift")+1);
1289         thermal_lift_start(ev);
1290         }
1291 else if (ev.lift_flag == 3)
1292         {
1293         ev.lift_r = cNode.getNode("thermal-lift").getValue();
1294         ev.n_entry_lift = getprop(lw~"effect-volumes/number-active-lift");      
1295         setprop(lw~"effect-volumes/number-active-lift",getprop(lw~"effect-volumes/number-active-lift")+1);
1296         wave_lift_start(ev);
1297         }
1298
1299 }
1300
1301
1302
1303 var effect_volume_stop = func (ev) {
1304
1305 var cNode = props.globals.getNode(lw~"current");
1306
1307
1308 if (ev.vis_flag == 1)
1309         {
1310
1311         var n_active = getprop(lw~"effect-volumes/number-active-vis");
1312
1313         
1314         var n_entry = ev.n_entry_vis;   
1315
1316         # if no other nodes affecting property are active, restore to outside
1317         # else restore settings as they have been when entering the volume when the number
1318         # of active volumes is the same as on entry (i.e. volumes are nested), otherwise
1319         # leave property at current because new definitions are already active and should not
1320         # be cancelled
1321         
1322         if (n_active ==1){var vis = interpolated_conditions.visibility_m;}
1323         else if ((n_active -1) == n_entry) 
1324                 {var vis = ev.vis_r;}
1325         else {var vis = cNode.getNode("visibility-m").getValue();}
1326         cNode.getNode("visibility-m").setValue(vis);
1327         compat_layer.setVisibilitySmoothly(vis);
1328         
1329         # and subtract from the counter
1330         setprop(lw~"effect-volumes/number-active-vis",getprop(lw~"effect-volumes/number-active-vis")-1);
1331         }
1332 if (ev.rain_flag == 1)
1333         {
1334         var n_active = getprop(lw~"effect-volumes/number-active-rain");
1335         var n_entry = ev.n_entry_rain;
1336
1337         if (n_active ==1){var rain = interpolated_conditions.rain_norm;}
1338         else if ((n_active -1) == n_entry)
1339                  {var rain = ev.rain_r;}
1340         else {var rain = cNode.getNode("rain-norm").getValue();}
1341         cNode.getNode("rain-norm").setValue(rain);
1342         compat_layer.setRain(rain);
1343         setprop(lw~"effect-volumes/number-active-rain",getprop(lw~"effect-volumes/number-active-rain")-1);
1344         }
1345
1346 if (ev.snow_flag == 1)
1347         {
1348         var n_active = getprop(lw~"effect-volumes/number-active-snow");
1349         var n_entry = ev.n_entry_snow;  
1350
1351         if (n_active ==1){var snow = interpolated_conditions.snow_norm;}
1352         else if ((n_active -1) == n_entry)
1353                 {var snow = ev.snow_r;}
1354         else {var snow = cNode.getNode("snow-norm").getValue();}
1355         cNode.getNode("snow-norm").setValue(snow);
1356         compat_layer.setSnow(snow);
1357         setprop(lw~"effect-volumes/number-active-snow",getprop(lw~"effect-volumes/number-active-snow")-1);
1358         }
1359
1360 if (ev.turb_flag == 1)
1361         {
1362         var n_active = getprop(lw~"effect-volumes/number-active-turb");
1363         var n_entry = ev.n_entry_turb;
1364         if (n_active ==1){var turbulence = interpolated_conditions.turbulence;}
1365         else if ((n_active -1) == n_entry) 
1366                  {var turbulence = ev.turb_r;}
1367         else {var turbulence = cNode.getNode("turbulence").getValue();}
1368         cNode.getNode("turbulence").setValue(turbulence);
1369         compat_layer.setTurbulence(turbulence);
1370         setprop(lw~"effect-volumes/number-active-turb",getprop(lw~"effect-volumes/number-active-turb")-1);
1371         }
1372
1373 if (ev.sat_flag == 1)
1374         {
1375         var n_active = getprop(lw~"effect-volumes/number-active-sat");
1376         var n_entry = ev.n_entry_sat;
1377         if (n_active ==1){var saturation = 1.0;}
1378         else if ((n_active -1) == n_entry) 
1379                  {var saturation = ev.sat_r;}
1380         else {var saturation = getprop("/rendering/scene/saturation");}
1381         compat_layer.setLightSmoothly(saturation);
1382         setprop(lw~"effect-volumes/number-active-sat",getprop(lw~"effect-volumes/number-active-sat")-1);
1383         }
1384
1385 if (ev.lift_flag == 1)
1386         {
1387         var n_active = getprop(lw~"effect-volumes/number-active-lift");
1388         var n_entry = ev.n_entry_lift;
1389         if (n_active ==1){var lift = interpolated_conditions.thermal_lift;}
1390         else if ((n_active -1) == n_entry)
1391                  {var lift = ev.lift_r;}
1392         else {var lift = cNode.getNode("thermal-lift").getValue();}
1393         cNode.getNode("thermal-lift").setValue(lift);
1394         compat_layer.setLift(lift);
1395         setprop(lw~"effect-volumes/number-active-lift",getprop(lw~"effect-volumes/number-active-lift")-1);
1396         }
1397 else if (ev.lift_flag == 2)
1398         {
1399         thermal_lift_stop();
1400         setprop(lw~"effect-volumes/number-active-lift",getprop(lw~"effect-volumes/number-active-lift")-1);
1401         }
1402 else if (ev.lift_flag == 3)
1403         {
1404         wave_lift_stop();
1405         setprop(lw~"effect-volumes/number-active-lift",getprop(lw~"effect-volumes/number-active-lift")-1);
1406         }
1407
1408 }
1409
1410
1411
1412 #########################################
1413 # compute thermal lift in detailed model 
1414 #########################################
1415
1416 var ts_factor = func (t, alt, height) {
1417
1418 var t1 = 0.1; # fractional time at which lift is fully developed 
1419 var t2 = 0.9; # fractional time at which lift starts to decay
1420 var t3 = 1.0; # fractional time at which lift is gone
1421
1422 # no time dependence modelled yet
1423 # return 1.0; 
1424
1425
1426
1427 var t_a = t - (alt/height) * t1 - t1;
1428
1429 if (t_a<0) {return 0.0;}
1430 else if (t_a<t1) {return 0.5 + 0.5 * math.cos((1.0-t_a/t1)* math.pi);}
1431 else if (t_a < t2) {return 1.0;}
1432 else {return 0.5 - 0.5 * math.cos((1.0-(t2-t_a)/(t3-t2))*math.pi);}
1433 }
1434
1435 var tl_factor = func (t, alt, height) {
1436
1437 var t1 = 0.1; # fractional time at which lift is fully developed 
1438 var t2 = 0.9; # fractional time at which lift starts to decay
1439 var t3 = 1.0; # fractional time at which lift is gone
1440
1441 # no time dependence modelled yet
1442 # return 1.0; 
1443
1444 var t_a = t - (alt/height) * t1;
1445
1446 if (t_a<0) {return 0.0;}
1447 else if (t_a<t1) {return 0.5 + 0.5 * math.cos((1.0-t_a/t1)* math.pi);}
1448 else if (t_a < t2) {return 1.0;}
1449 else  {return 0.5 - 0.5 * math.cos((1.0-(t2-t_a)/(t3-t2))*math.pi);}      
1450 }
1451
1452
1453 var calcLift_max = func (alt, max_lift, height) {
1454     
1455 alt_agl = getprop("/position/altitude-agl-ft");
1456
1457 # no lift below ground
1458 if (alt_agl < 0.0) {return 0.0;}   
1459     
1460 # lift ramps up to full within 200 m
1461 else if (alt_agl < 200.0*m_to_ft) 
1462         {return max_lift * 0.5 * (1.0 + math.cos((1.0-alt_agl/(200.0*m_to_ft))*math.pi));}
1463
1464 # constant max. lift in main body
1465 else if ((alt_agl > 200.0*m_to_ft) and (alt < height))
1466         {return max_lift;}
1467
1468 # decreasing lift from cloudbase to 10% above base
1469 else if ((alt > height ) and (alt < height*1.1)) 
1470         {return max_lift * 0.5 * (1.0 - math.cos((1.0-10.0*(alt-height)/height)*math.pi));}
1471     
1472 # no lift available above
1473 else {return 0.0;}
1474 }
1475
1476
1477
1478 var calcLift = func (d, alt, R, height, cn, sh, max_lift, f_lift_radius, t) {
1479
1480 # radius of slice at given altitude
1481 var r_total = (cn + alt/height*(1.0-cn)) * (R - R * (1.0- sh ) * (1.0 - ((2.0*alt/height)-1.0)*((2.0*alt/height)-1.0)));
1482
1483
1484 # print("r_total: ", r_total, "d: ",d);
1485 # print("alt: ", alt, "height: ",height);
1486
1487 # no lift if we're outside the radius or above the thermal
1488 if ((d > r_total) or (alt > 1.1*height)) { return 0.0; } 
1489
1490 # fraction of radius providing lift
1491 var r_lift = f_lift_radius * r_total;
1492
1493 # print("r_lift: ", r_lift);
1494
1495 # if we are in the sink portion, get the max. sink for this time and altitude and adjust for actual position
1496 if ((d < r_total ) and (d > r_lift)) 
1497         {
1498         var s_max = 0.5 * calcLift_max(alt, max_lift, height) * ts_factor(t, alt, height);
1499         # print("s_max: ", s_max);
1500         return s_max * math.sin(math.pi * (1.0 + (d-r_lift) * (1.0/(r_total - r_lift))));
1501         }
1502 # else we are in the lift portion, get the max. lift for this time and altitude and adjust for actual position
1503 else
1504         {  
1505         var l_max = calcLift_max(alt, max_lift, height) * tl_factor(t, alt, height);
1506         # print("l_max: ", l_max);
1507         return l_max * math.cos(math.pi * (d/(2.0 * r_lift)));
1508         }
1509 }
1510
1511 #########################################
1512 # compute wave lift in detailed model 
1513 #########################################
1514
1515 var calcWaveLift = func (x,y, alt) {
1516
1517 var lift = wave.max_lift * math.cos((y/wave.y) * 1.5 * math.pi);
1518
1519 if (abs(x)/wave.x > 0.9)
1520         {
1521         lift = lift * (abs(x) - 0.9 * wave.x)/(0.1 * wave.x); 
1522         }
1523
1524
1525
1526 lift = lift * 2.71828 * math.exp(-alt/wave.height) * alt/wave.height;
1527
1528 var alt_agl = getprop("/position/altitude-agl-ft");
1529
1530 if (alt_agl < 1000.0)
1531         {
1532         lift = lift * (alt_agl/1000.0) * (alt_agl/1000.0);
1533         }
1534
1535 return lift;
1536 }
1537         
1538
1539
1540
1541
1542
1543 ###########################################################
1544 # place a single cloud into a vector to be processed
1545 # separately
1546 ###########################################################
1547
1548 var create_cloud_vec = func(path, lat, lon, alt, heading) {
1549
1550 if (path == "new") # we have to switch to new cloud generating routines
1551         {
1552         local_weather.cloudAssembly.lat = lat;
1553         local_weather.cloudAssembly.lon = lon;
1554         local_weather.cloudAssembly.alt = alt;  
1555         local_weather.cloudAssembly.top_shade = top_shade;
1556
1557         #print(lat," ",long, " ", alt);
1558
1559         if (dynamics_flag == 1)
1560                 {
1561                 local_weather.cloudAssembly.mean_alt = cloud_mean_altitude;
1562                 local_weather.cloudAssembly.flt = cloud_fractional_lifetime;
1563                 local_weather.cloudAssembly.evolution_timestamp = cloud_evolution_timestamp;
1564                 local_weather.cloudAssembly.rel_alt = cloudAssembly.alt - cloud_mean_altitude;
1565                 }
1566
1567         append(cloudAssemblyArray,cloudAssembly);
1568
1569         # at this point we insert tracers for the depth buffer
1570         
1571         #if (local_weather.cloudAssembly.tracer_flag == 1)
1572         #       {       
1573         #       tracerAssembly = local_weather.cloud.new("Tracer", "default");
1574         #       tracerAssembly.texture_sheet = "/Models/Weather/nimbus_sheet1.rgb";
1575         #       tracerAssembly.n_sprites = 1;
1576         #       tracerAssembly.bottom_shade = 0.0;
1577         #       tracerAssembly.top_shade = 0.0;
1578         #       tracerAssembly.num_tex_x = 1;
1579         #       tracerAssembly.num_tex_y = 1;
1580         #       tracerAssembly.lat = lat;
1581         #       tracerAssembly.lon = lon;
1582         #       tracerAssembly.alt = alt + local_weather.cloudAssembly.min_height *0.35 * m_to_ft ;
1583         #       tracerAssembly.min_width = local_weather.cloudAssembly.min_width * 0.35;
1584         #       tracerAssembly.max_width = local_weather.cloudAssembly.max_width * 0.35;
1585         #       tracerAssembly.min_height = local_weather.cloudAssembly.min_height * 0.35;
1586         #       tracerAssembly.max_height = local_weather.cloudAssembly.max_height * 0.35;
1587         #       tracerAssembly.min_cloud_width = local_weather.cloudAssembly.min_cloud_width * 0.35;
1588         #       tracerAssembly.min_cloud_height = local_weather.cloudAssembly.min_cloud_height * 0.35;
1589         #       tracerAssembly.z_scale = local_weather.cloudAssembly.z_scale;
1590         #       append(cloudAssemblyArray,tracerAssembly);
1591         #       }
1592
1593         return;
1594         }
1595
1596 append(clouds_path,path);
1597 append(clouds_lat,lat);
1598 append(clouds_lon,lon);
1599 append(clouds_alt,alt);
1600 append(clouds_orientation,heading);
1601
1602 # globals (needed for Cumulus clouds) should be set if needed by the main cloud generating call
1603
1604 if (dynamics_flag ==1)
1605         {
1606         append(clouds_mean_alt, cloud_mean_altitude);
1607         append(clouds_flt, cloud_fractional_lifetime);
1608         append(clouds_evolution_timestamp,cloud_evolution_timestamp);
1609         }
1610
1611 }
1612 ###########################################################
1613 # clear all clouds and effects
1614 ###########################################################
1615
1616 var clear_all = func {
1617
1618 # clear the clouds and models
1619
1620 var cloudNode = props.globals.getNode(lw~"clouds", 1);
1621 cloudNode.removeChildren("tile");
1622
1623 var modelNode = props.globals.getNode("models", 1).getChildren("model");
1624
1625 foreach (var m; modelNode)
1626         {
1627         var l = m.getNode("tile-index",1).getValue();
1628         if (l != nil)
1629                 {
1630                 m.remove();
1631                 }
1632         }
1633
1634
1635 # remove the hard-coded clouds
1636
1637 foreach (c; weather_tile_management.cloudArray)
1638         {
1639         c.remove();
1640         }
1641 setsize(weather_tile_management.cloudArray,0);
1642
1643 # reset pressure continuity
1644
1645 weather_tiles.last_pressure = 0.0;
1646
1647 # stop all loops
1648
1649 setprop(lw~"effect-loop-flag",0);
1650 setprop(lw~"interpolation-loop-flag",0);
1651 setprop(lw~"tile-loop-flag",0);
1652 setprop(lw~"lift-loop-flag",0);
1653 setprop(lw~"wave-loop-flag",0);
1654 setprop(lw~"dynamics-loop-flag",0);
1655 setprop(lw~"timing-loop-flag",0);
1656 setprop(lw~"buffer-loop-flag",0);
1657 setprop(lw~"housekeeping-loop-flag",0);
1658 setprop(lw~"convective-loop-flag",0);
1659
1660 weather_dynamics.convective_loop_kill_flag = 1; # long-running loop needs a different scheme to end
1661
1662 # also remove rain, snow, haze and light effects
1663
1664 compat_layer.setRain(0.0);
1665 compat_layer.setSnow(0.0);
1666 compat_layer.setLight(1.0);
1667
1668
1669 # set placement indices to zero
1670
1671 setprop(lw~"clouds/placement-index",0);
1672 setprop(lw~"clouds/model-placement-index",0);
1673 setprop(lw~"effect-volumes/effect-placement-index",0);
1674 setprop(lw~"effect-volumes/number",0);
1675 setprop(lw~"effect-volumes/number-active-rain",0);
1676 setprop(lw~"effect-volumes/number-active-snow",0);
1677 setprop(lw~"effect-volumes/number-active-vis",0);
1678 setprop(lw~"effect-volumes/number-active-turb",0);
1679 setprop(lw~"effect-volumes/number-active-lift",0);
1680 setprop(lw~"effect-volumes/number-active-sat",0);
1681 setprop(lw~"tiles/tile-counter",0);
1682
1683
1684 # remove any quadtrees and arrays
1685
1686 settimer ( func { setsize(weather_dynamics.cloudQuadtrees,0);},0.1); # to avoid error generation in this frame
1687 setsize(effectVolumeArray,0);
1688 n_effectVolumeArray = 0;
1689
1690
1691 # clear any wxradar echos
1692
1693 if (wxradar_support_flag ==1)
1694         {props.globals.getNode("/instrumentation/wxradar", 1).removeChildren("storm");}
1695
1696 # if we have used METAR, we may no longer want to do so
1697
1698 metar_flag = 0;
1699
1700
1701 settimer ( func {
1702         setsize(weather_tile_management.modelArrays,0);
1703         setsize(weather_dynamics.tile_wind_direction,0);
1704         setsize(weather_dynamics.tile_wind_speed,0);
1705         setsize(weather_tile_management.cloudBufferArray,0);
1706         setsize(weather_tile_management.cloudSceneryArray,0);
1707         setsize(alt_20_array,0);
1708         setsize(alt_50_array,0);
1709         setsize(alt_min_array,0);
1710         setsize(alt_mean_array,0);
1711         setsize(weather_dynamics.tile_convective_altitude,0);
1712         setsize(weather_dynamics.tile_convective_strength,0);
1713         setsize(weatherStationArray,0);
1714         setsize(windIpointArray,0);
1715         setsize(atmosphereIpointArray,0);
1716         setprop(lw~"clouds/buffer-count",0);
1717         setprop(lw~"clouds/cloud-scenery-count",0);
1718         weather_tile_management.n_cloudSceneryArray = 0;
1719         compat_layer.setScattering(0.8);
1720         compat_layer.setOvercast(0.0);
1721         setprop(lwi~"ipoint-number",0);
1722         setprop(lwi~"atmosphere-ipoint-number", 0);
1723         },0);
1724
1725 setprop(lw~"tmp/presampling-status", "idle");
1726
1727 # reset the random store
1728
1729 weather_tiles.rnd_store = rand();
1730
1731 # default 3d clouds layer wrapping back on, just in case
1732
1733 setprop("/sim/rendering/clouds3d-wrap",1);
1734
1735 # indicate that we are no longer running
1736
1737
1738 local_weather_running_flag = 0;
1739
1740 }
1741
1742
1743
1744 ###########################################################
1745 # detailed Cumulus clouds created from multiple cloudlets
1746 ###########################################################
1747
1748 var create_detailed_cumulus_cloud = func (lat, lon, alt, size) {
1749
1750
1751 # various distribution biases
1752
1753 var edge_bias = convective_texture_mix;
1754 size = size + convective_size_bias;
1755 height_bias = 1.0;
1756 if (edge_bias > 0.0) {height_bias = height_bias +  15.0 *edge_bias + 20.0 * rand() * edge_bias;}
1757
1758
1759 #height_bias = 6.0;
1760
1761         
1762         if (size > 2.0)
1763                 {
1764                 if (rand() > (size - 2.0))
1765                         {create_cumulonimbus_cloud(lat, lon, alt, size); }
1766                 else
1767                         {create_cumulonimbus_cloud_rain(lat, lon, alt, size, 0.1 + 0.2* rand());}
1768                 return;
1769                 }
1770
1771         else if (size>1.5)
1772                 {
1773                 var type = "Congestus";
1774
1775                 var height = 400;
1776                 var n = 3;
1777                 var x = 700.0;
1778                 var y = 200.0;
1779                 var edge = 0.2;
1780                 
1781                 var alpha = rand() * 180.0;
1782                 edge = edge + edge_bias;                
1783
1784                 create_streak(type,lat,lon, alt+ 0.3* (height )-offset_map["Congestus"], height,n,0.0,edge,x,1,0.0,0.0,y,alpha,1.0);
1785
1786                 var type = "Cu (volume)";
1787                 var height = 400;
1788                 var n = 10 + int(height_bias);
1789                 var x = 1400.0;
1790                 var y = 400.0;
1791                 var edge = 0.2;
1792                 
1793                 edge = edge + edge_bias;                
1794
1795                 create_streak(type,lat,lon, alt+ 0.5* (height * height_bias )-offset_map["Cumulus"], height * height_bias ,n,0.0,edge,x,1,0.0,0.0,y,alpha,1.0);
1796
1797                 var btype = "Congestus bottom";
1798                 var n_b = 6;
1799                 height_bias = 1.0;
1800
1801                 create_streak(btype,lat,lon, alt -offset_map["Congestus"] -900.0, 100.0,n_b,0.0,edge,0.3*x,1,0.0,0.0,0.3*y,alpha,1.0);
1802
1803                 }
1804         else if (size>1.1)
1805                 {
1806                 var type = "Cumulus (cloudlet)";
1807                 var btype = "Cumulus bottom";
1808                 var height = 200;
1809                 var n = 6 + int(height_bias);
1810                 var n_b = 2;
1811                 var x = 900.0;
1812                 var y = 200.0;
1813                 var edge = 0.2;
1814
1815                 var alpha = rand() * 180.0;
1816                 edge = edge + edge_bias;
1817                 create_streak(type,lat,lon, alt+ 0.5* (height* height_bias )-offset_map["Cumulus"], height * height_bias,n,0.0,edge,x,1,0.0,0.0,y,alpha,1.0);
1818
1819                 height_bias = 1.0;
1820                 create_streak(btype,lat,lon, alt -offset_map["Cumulus"] - 200.0, 100.0,n_b,0.0,edge,0.3*x,1,0.0,0.0,0.3*y,alpha,1.0);
1821
1822                 }
1823         else if (size>0.8)
1824                 {
1825                 var type = "Cumulus (cloudlet)";
1826                 var height = 150;
1827                 var n = 4 + int(height_bias);
1828                 var x = 300.0;
1829                 var y = 300.0;
1830                 var edge = 0.3;
1831
1832                 var alpha = rand() * 180.0;
1833                 edge = edge + edge_bias;
1834                 create_streak(type,lat,lon, alt+ 0.5* (height * height_bias )-offset_map["Cumulus"], height * height_bias,n,0.0,edge,x,1,0.0,0.0,y,alpha,1.0);
1835
1836                 n = 2;
1837                 x = 700.0;
1838                 y = 200.0;
1839                 edge = 1.0;
1840                 create_streak(type,lat,lon, alt+ 0.5* (height*height_bias )-offset_map["Cumulus"], height * height_bias,n,0.0,edge,x,1,0.0,0.0,y,alpha,1.0);
1841                 }
1842
1843         else if (size>0.4)
1844                 {
1845                 var type = "Cumulus (cloudlet)";
1846                 var height = 100;
1847                 var n = 2 + int(height_bias * 0.5);
1848                 var x = 600.0;
1849                 var y = 100.0;
1850                 var edge = 1.0;
1851
1852                 var alpha = rand() * 180.0;
1853                 edge = edge + edge_bias;
1854                 create_streak(type,lat,lon, alt+ 0.5* (height * height_bias)-offset_map["Cumulus"], height * height_bias,n,0.0,edge,x,1,0.0,0.0,y,alpha,1.0);
1855                 }
1856         else 
1857                 {
1858                 var type = "Cumulus (whisp)";
1859                 var height = 100;
1860                 var n = 1;
1861                 var x = 100.0;
1862                 var y = 100.0;
1863                 var edge = 1.0;
1864
1865                 var alpha = rand() * 180.0;
1866                 edge = edge + edge_bias;
1867                 create_streak(type,lat,lon, alt+ 0.3* (height )-offset_map["Cumulus"], height,n,0.0,edge,x,1,0.0,0.0,y,alpha,1.0);
1868                 }
1869
1870
1871 ###########################################################
1872 # detailed small Cumulonimbus clouds created from multiple cloudlets
1873 ###########################################################
1874
1875 var create_cumulonimbus_cloud = func(lat, lon, alt, size) {
1876
1877
1878 create_cloudbox("Cb_box", lat, lon, alt, 2500.0,2000.0, 1000.0,10, 0.2, 0.1, 1.0, 1, 0.8, 0.1, 6);
1879
1880 #create_cloudbox = func (type, blat, blon, balt, dx,dy,dz,n, f_core, r_core, h_core, n_core, f_bottom, h_bottom, n_bottom)
1881 }
1882
1883 ###########################################################
1884 # detailed small Cumulonimbus and rain created from multiple cloudlets
1885 ###########################################################
1886
1887 var create_cumulonimbus_cloud_rain = func(lat, lon, alt, size, rain) {
1888
1889
1890 create_cloudbox("Cb_box", lat, lon, alt, 2500.0,2000.0, 1000.0,10, 0.2, 0.1, 1.0, 1, 0.8, 0.1, 6);
1891
1892
1893
1894 # place a rain texture
1895
1896 var path = "Models/Weather/rain2.xml";
1897 if (thread_flag == 1)
1898                                 {create_cloud_vec(path, lat, lon, alt, 0.0);}
1899                         else 
1900                                 {compat_layer.create_cloud(path, lat, lon, alt, 0.0);}
1901
1902                 
1903
1904 # and some rain underneath
1905
1906 create_effect_volume(1, lat, lon, 2000.0, 2000.0, 0.0, 0.0, alt+1000.0, 8000.0 + 8000.0 * rand(), rain, -1, -1, -1 ,1,-1 );
1907
1908
1909 }
1910
1911
1912 ###########################################################
1913 # wrappers for convective cloud system to distribute
1914 # call across several frames if needed
1915 ###########################################################
1916
1917 var create_cumosys = func (blat, blon, balt, nc, size) {
1918
1919
1920 # realistic Cumulus has somewhat larger models, so compensate to get the same coverage
1921 if (detailed_clouds_flag == 1) 
1922         {nc = int(0.7 * nc);}
1923
1924 nc = int(nc / cumulus_efficiency_factor);
1925
1926 if (thread_flag ==  1)
1927         {setprop(lw~"tmp/convective-status", "computing");
1928         cumulus_loop(blat, blon, balt, nc, size);}
1929
1930 else
1931         {create_cumulus(blat, blon, balt, nc, size);
1932         if (debug_output_flag == 1) 
1933                 {print("Convective system done!");}
1934         }
1935 }
1936
1937
1938
1939 var cumulus_loop = func (blat, blon, balt, nc, size) {
1940
1941 if (local_weather_running_flag == 0) {return;}
1942
1943 if (local_weather.features.fast_geodinfo == 0)
1944         {var n = int(25/cumulus_efficiency_factor);}
1945 else
1946         {var n = int(200/cumulus_efficiency_factor);}
1947
1948 if (nc < 0) 
1949         {
1950         if (debug_output_flag == 1) 
1951                 {print("Convective system done!");}
1952         setprop(lw~"tmp/convective-status", "idle");
1953         assemble_effect_array();
1954         convective_size_bias = 0.0;
1955         height_bias = 1.0;
1956         return;
1957         }
1958
1959 create_cumulus(blat, blon, balt, n, size);
1960
1961 settimer( func {cumulus_loop(blat, blon, balt, nc-n, size) },0);
1962 }
1963
1964 ###########################################################
1965 # place a convective cloud system
1966 ###########################################################
1967
1968 var create_cumulus = func (blat, blon, balt, nc, size) {
1969
1970
1971
1972 var path = "Models/Weather/blank.ac";
1973 var i = 0;
1974 var p = 0.0;
1975 var rn = 0.0;
1976 var place_lift_flag = 0;
1977 var strength = 0.0;
1978 var detail_flag = detailed_clouds_flag;
1979
1980 var alpha = getprop(lw~"tmp/tile-orientation-deg") * math.pi/180.0; # the tile orientation
1981
1982 var tile_index = getprop(lw~"tiles/tile-counter");
1983 var alt_base = alt_20_array[tile_index -1];
1984
1985 #if (detailed_terrain_interaction_flag == 1)
1986 #       {
1987         #var tile_index = getprop(lw~"tiles/tile-counter");
1988         #var alt_min = alt_min_array[tile_index-1];
1989         #var alt_mean = alt_mean_array[tile_index -1];
1990         #var alt_median = alt_50_array[tile_index -1];
1991         #var alt_base = alt_20_array[tile_index -1];
1992         #var alt_min = getprop(lw~"tmp/tile-alt-min-ft");
1993         #var alt_mean = getprop(lw~"tmp/tile-alt-mean-ft");
1994         #var alt_median = getprop(lw~"tmp/tile-alt-median-ft");
1995         #var alt_base = getprop(lw~"tmp/tile-alt-offset-ft");
1996 #       }
1997
1998 var sec_to_rad = 2.0 * math.pi/86400; # conversion factor for sinusoidal dependence on daytime
1999
2000 calc_geo(blat);
2001
2002 # get the local time of the day in seconds
2003
2004 var t = getprop("sim/time/utc/day-seconds");
2005 t = t + getprop("sim/time/local-offset");
2006
2007 # print("t is now:", t);
2008
2009 # and make a simple sinusoidal model of thermal strength
2010
2011 # daily variation in number of thermals, peaks at noon
2012 var t_factor1 = 0.5 * (1.0-math.cos((t * sec_to_rad))); 
2013
2014 # daily variation in strength of thermals, peaks around 15:30
2015 var t_factor2 = 0.5 * (1.0-math.cos((t * sec_to_rad)-0.9)); 
2016
2017
2018 # number of possible thermals equals overall strength times daily variation times geographic variation
2019 # this is a proxy for solar thermal energy
2020
2021 nc = t_factor1 * nc * math.cos(blat/180.0*math.pi); 
2022
2023 # var thermal_conditions = getprop(lw~"config/thermal-properties");
2024
2025
2026 while (i < nc) {
2027
2028         p = 0.0;
2029         place_lift_flag = 0;
2030         strength = 0.0;
2031
2032         # pick a trial position inside the tile and rotate by tile orientation angle
2033         var x = (2.0 * rand() - 1.0) * size;
2034         var y = (2.0 * rand() - 1.0) * size; 
2035
2036         var lat = blat + (y * math.cos(alpha) - x * math.sin(alpha)) * m_to_lat;
2037         var lon = blon + (x * math.cos(alpha) + y * math.sin(alpha)) * m_to_lon;
2038
2039         # now check ground cover type on chosen spot
2040         var info = geodinfo(lat, lon);
2041
2042         if (info != nil) {
2043         var elevation = info[0] * m_to_ft;
2044         if (info[1] != nil){
2045          var landcover = info[1].names[0];
2046          if (contains(landcover_map,landcover)) {p = p + landcover_map[landcover];}
2047          else {print(p, " ", info[1].names[0]);}
2048         }}
2049         else {
2050                 # to avoid gaps, we create default clouds
2051
2052                 p = p + 0.1;            
2053                 var elevation = alt_base;
2054                 # continue;
2055                 }
2056
2057
2058         # apply some optional corrections, biases clouds towards higher elevations
2059
2060         var terrain_altitude_factor = 1.0;
2061         var terrain_strength_factor = 1.0;
2062
2063         if (detailed_terrain_interaction_flag == 1)
2064                 {
2065                 
2066                 terrain_altitude_factor = get_terrain_altitude_factor(tile_index, balt, elevation);
2067                 terrain_strength_factor = get_terrain_strength_factor(terrain_altitude_factor);
2068
2069                 }
2070
2071
2072         # then decide if the thermal energy at the spot generates an updraft and a cloud
2073
2074         if (rand() < (p * cumulus_efficiency_factor * terrain_altitude_factor)) # we decide to place a cloud at this spot
2075                 {
2076         
2077
2078                 # check if we have a terrain elevation analysis available and can use a 
2079                 # detailed placement altitude correction
2080
2081                 if (presampling_flag == 1) 
2082                         {
2083                         
2084                         if (detailed_terrain_interaction_flag == 1)
2085                                 {
2086                                 var grad = get_terrain_gradient(lat, lon, elevation, alpha, 1000.0);
2087                                 }
2088                         else 
2089                                 {var grad = 0.0;}
2090
2091
2092                         var place_alt = get_convective_altitude(balt, elevation, getprop(lw~"tiles/tile-counter"), grad);
2093                         }
2094                 else {var place_alt = balt;}
2095                 
2096                 # no cloud placement into the ground
2097                 if (place_alt < elevation) {continue;}
2098
2099                 # if we're in a lee, we may not want to place the cloud
2100
2101                 if (detailed_terrain_interaction_flag == 1)
2102                                 {
2103                                 var p_lee_suppression = get_lee_bias(grad, tile_index);
2104                                 if (rand() > p_lee_suppression) {continue;} 
2105                                 }
2106
2107         
2108                 # now decide on the strength of the thermal at the spot and on cloud size
2109
2110                 var rn = rand();
2111                 strength = (1.5 * rn + (2.0 * p * terrain_strength_factor)) * t_factor2;  
2112                 
2113                 # the terrain effect cannot create Cb development, so we have to curb
2114                 # the strength if it would not have been Cb otherwise
2115
2116                 if (strength > 2.0)
2117                         {
2118                         if (((1.5 * rn + (2.0 * p)) * t_factor2) < 2.0)
2119                                 {strength = 1.7 + rand() * 0.2;}
2120                         }
2121
2122
2123                 if (strength > 1.0) {place_lift_flag = 1;}
2124
2125                 cloud_mean_altitude = place_alt;
2126                 cloud_fractional_lifetime = rand();
2127                 cloud_evolution_timestamp = weather_dynamics.time_lw;
2128
2129
2130
2131                 if (generate_thermal_lift_flag != 3) # no clouds if we produce blue thermals
2132                         {               
2133                         create_detailed_cumulus_cloud(lat, lon, place_alt, strength);
2134                         }       
2135
2136                 # now see if we need to create a thermal - first check the flag
2137                 if (generate_thermal_lift_flag == 1) # thermal by constant
2138                         {
2139                         # now check if convection is strong
2140                         if (place_lift_flag == 1)
2141                                 {
2142                                 var lift = 3.0 + 10.0 * (strength -1.0);
2143                                 var radius = 500 + 500 * rand();
2144                                 #print("Lift: ", lift * ft_to_m - 1.0);
2145                                 create_effect_volume(1, lat, lon, radius, radius, 0.0, 0.0, place_alt+500.0, -1, -1, -1, -1, lift, 1,-1);
2146                                 } # end if place_lift_flag              
2147                         } # end if generate-thermal-lift-flag   
2148                 else if ((generate_thermal_lift_flag == 2) or (generate_thermal_lift_flag == 3)) # thermal by function
2149                         {
2150
2151                         if (place_lift_flag == 1)
2152                                 {
2153                                 var lift = (3.0 + 10.0 * (strength -1.0))/thermal_conditions;
2154                                 var radius = (500 + 500 * rand())*thermal_conditions;
2155
2156                                 create_effect_volume(1, lat, lon, 1.1*radius, 1.1*radius, 0.0, 0.0, place_alt*1.15, -1, -1, -1, lift*0.03, lift, -2,-1);
2157                                 } # end if place_lift_flag
2158
2159                         } # end if generate-thermal-lift-flag
2160
2161
2162                 } # end if rand < p
2163         i = i + 1;
2164         } # end while
2165
2166 }
2167
2168
2169
2170
2171
2172 #################################################################
2173 # respawn convective clouds to compensate for decay
2174 # the difference being that new clouds get zero fractional 
2175 # lifetime and are placed based on terrain with a different weight
2176 ##################################################################
2177
2178 var recreate_cumulus = func (blat, blon, balt, alpha, nc, size, tile_index) {
2179
2180 var path = "Models/Weather/blank.ac";
2181 var i = 0;
2182 var p = 0.0;
2183 var rn = 0.0;
2184 var place_lift_flag = 0;
2185 var strength = 0.0;
2186 var detail_flag = detailed_clouds_flag;
2187
2188 alpha = alpha * math.pi/180.0; # the tile orientation
2189
2190 var sec_to_rad = 2.0 * math.pi/86400; # conversion factor for sinusoidal dependence on daytime
2191
2192 # current aircraft position
2193
2194 var alat = getprop("position/latitude-deg");
2195 var alon = getprop("position/longitude-deg");
2196
2197 # get the local time of the day in seconds
2198
2199 var t = getprop("sim/time/utc/day-seconds");
2200 t = t + getprop("sim/time/local-offset");
2201
2202
2203 # and make a simple sinusoidal model of thermal strength
2204
2205 # daily variation in number of thermals, peaks at noon
2206 var t_factor1 = 0.5 * (1.0-math.cos((t * sec_to_rad))); 
2207
2208 # daily variation in strength of thermals, peaks around 15:30
2209 var t_factor2 = 0.5 * (1.0-math.cos((t * sec_to_rad)-0.9)); 
2210
2211
2212 # number of possible thermals equals overall strength times daily variation times geographic variation
2213 # this is a proxy for solar thermal energy
2214
2215 nc = t_factor1 * nc * math.cos(blat/180.0*math.pi); 
2216
2217 # var thermal_conditions = getprop(lw~"config/thermal-properties");
2218
2219 var alt_base = alt_20_array[tile_index -1];
2220
2221 while (i < nc) {
2222
2223         p = 0.0;
2224         place_lift_flag = 0;
2225         strength = 0.0;
2226
2227         # pick a trial position inside the tile and rotate by tile orientation angle
2228         var x = (2.0 * rand() - 1.0) * size;
2229         var y = (2.0 * rand() - 1.0) * size; 
2230
2231         var lat = blat + (y * math.cos(alpha) - x * math.sin(alpha)) * m_to_lat;
2232         var lon = blon + (x * math.cos(alpha) + y * math.sin(alpha)) * m_to_lon;
2233
2234         # check if the cloud would be spawned in visual range, if not don't bother
2235         var d_sq = calc_d_sq(alat, alon, lat, lon);
2236
2237         if (math.sqrt(d_sq) > weather_tile_management.cloud_view_distance)
2238                 {i = i+1; continue;}
2239
2240         # now check ground cover type on chosen spot
2241         var info = geodinfo(lat, lon);
2242
2243         if (info != nil) {
2244         var elevation = info[0] * m_to_ft;
2245         if (info[1] != nil){
2246          var landcover = info[1].names[0];
2247          if (contains(landcover_map,landcover)) {p = p + landcover_map[landcover];}
2248          else {print(p, " ", info[1].names[0]);}
2249         }}
2250         else {
2251                 # to avoid gaps, we create default clouds
2252
2253                 p = p + 0.1;            
2254                 var elevation = alt_base;
2255                 # continue;
2256                 }
2257
2258
2259         # apply some optional corrections, biases clouds towards higher elevations
2260
2261         var terrain_altitude_factor = 1.0;
2262         var terrain_strength_factor = 1.0;
2263
2264         if (detailed_terrain_interaction_flag == 1)
2265                 {
2266                 terrain_altitude_factor = get_terrain_altitude_factor(tile_index, balt, elevation);
2267                 terrain_strength_factor = get_terrain_strength_factor(terrain_altitude_factor);
2268                 }
2269
2270
2271
2272
2273         # check if to place a cloud with weight sqrt(p), the lifetime gets another sqrt(p) factor
2274         
2275         if (rand() > math.sqrt(p * cumulus_efficiency_factor * terrain_altitude_factor))
2276                 {i=i+1; continue;}
2277
2278
2279         # then calculate the strength of the updraft
2280                 
2281         strength = (1.5 * rand() + (2.0 * p * terrain_strength_factor)) * t_factor2; # the strength of thermal activity at the spot
2282         if (strength > 1.0)  
2283                 {
2284                 path = select_cloud_model("Cumulus","large"); place_lift_flag = 1;
2285                 }
2286         else {path = select_cloud_model("Cumulus","small");}
2287
2288         if (presampling_flag == 1) 
2289                 {
2290                 var place_alt = get_convective_altitude(balt, elevation, tile_index,0.0);
2291                 }
2292         else {var place_alt = balt;}
2293
2294
2295         # no cloud placement into the ground
2296         if (place_alt < elevation) {continue;}
2297
2298         # if we're in a lee, we may not want to place the cloud
2299
2300         if (detailed_terrain_interaction_flag == 1)
2301                         {
2302                         var p_lee_suppression = get_lee_bias(grad, tile_index);
2303                                 if (rand() > math.sqrt(p_lee_suppression)) {continue;} 
2304                         }
2305                 
2306         cloud_mean_altitude = place_alt;
2307         cloud_fractional_lifetime = 0.0;
2308         cloud_evolution_timestamp = weather_dynamics.time_lw;
2309
2310         compat_layer.cloud_mean_altitude = place_alt;
2311         compat_layer.cloud_flt = cloud_fractional_lifetime;
2312         compat_layer.cloud_evolution_timestamp = cloud_evolution_timestamp;
2313
2314         if (generate_thermal_lift_flag != 3) # no clouds if we produce blue thermals
2315                 {               
2316                 if (thread_flag == 1)
2317                         {
2318                         thread_flag = 0; # create clouds immediately
2319                         if (detail_flag == 0){compat_layer.create_cloud(path,lat,lon, place_alt, 0.0);}
2320                         else {create_detailed_cumulus_cloud(lat, lon, place_alt, strength);}
2321                         thread_flag = 1; # and restore threading
2322                         }
2323                 else
2324                         {
2325                         if (detail_flag == 0){compat_layer.create_cloud(path, lat, lon, place_alt, 0.0);}
2326                         else {create_detailed_cumulus_cloud(lat, lon, place_alt, strength);}
2327                         }
2328                 }       
2329
2330         if (generate_thermal_lift_flag == 1) # thermal by constant
2331                 {
2332                 if (place_lift_flag == 1)
2333                         {
2334                         var lift = 3.0 + 10.0 * (strength -1.0);
2335                         var radius = 500 + 500 * rand();
2336                         create_effect_volume(1, lat, lon, radius, radius, 0.0, 0.0, place_alt+500.0, -1, -1, -1, -1, lift, 1,-1);
2337                         } # end if place_lift_flag              
2338                 } # end if generate-thermal-lift-flag   
2339         else if ((generate_thermal_lift_flag == 2) or (generate_thermal_lift_flag == 3)) # thermal by function
2340                 {
2341                         if (place_lift_flag == 1)
2342                         {
2343                         var lift = (3.0 + 10.0 * (strength -1.0))/thermal_conditions;
2344                         var radius = (500 + 500 * rand())*thermal_conditions;
2345
2346                         create_effect_volume(1, lat, lon, 1.1*radius, 1.1*radius, 0.0, 0.0, place_alt*1.15, -1, -1, -1, lift*0.03, lift, -2,-1);
2347                         } # end if place_lift_flag
2348
2349                 } # end if generate-thermal-lift-flag
2350
2351
2352         i = i + 1;
2353         } # end while
2354
2355 }
2356
2357
2358
2359
2360
2361
2362
2363 ###########################################################
2364 # place a barrier cloud system 
2365 ###########################################################
2366
2367 var create_rise_clouds = func (blat, blon, balt, nc, size, winddir, dist) {
2368
2369 var path = "Models/Weather/blank.ac";
2370 var i = 0;
2371 var p = 0.0;
2372 var rn = 0.0;
2373 var nsample = 10;
2374 var counter = 0;
2375 var dir = (winddir + 180.0) * math.pi/180.0;
2376 var step = dist/nsample;
2377
2378 calc_geo(blat);
2379
2380 while (i < nc) {
2381
2382         counter = counter + 1;
2383         p = 0.0; 
2384
2385         var x = (2.0 * rand() - 1.0) * size;
2386         var y = (2.0 * rand() - 1.0) * size; 
2387
2388         var lat = blat + y * m_to_lat;
2389         var lon = blon + x * m_to_lon;
2390
2391         var elevation = compat_layer.get_elevation(lat, lon);
2392         
2393         #print("elevation: ", elevation, "balt: ", balt);
2394
2395         if ((elevation < balt) and (elevation != -1.0))
2396         {
2397         for (var j = 0; j<nsample; j=j+1)
2398                 {
2399                 d = j * step;
2400                 x = d * math.sin(dir);
2401                 y = d * math.cos(dir);
2402                 var tlat = lat + y * m_to_lat;
2403                 var tlon = lon + x * m_to_lon;
2404                 
2405                 #print("x: ", x, "y: ", y);
2406
2407                 var elevation1 = compat_layer.get_elevation(tlat,tlon); 
2408                 #print("elevation1: ", elevation1, "balt: ", balt);
2409         
2410                 if (elevation1 > balt)
2411                         {
2412                         p = 1.0 - j * (1.0/nsample);
2413                         #p = 1.0;
2414                         break;
2415                         }
2416                 
2417                 }
2418         }
2419         if (counter > 500) {print("Cannot place clouds - exiting..."); i = nc;}
2420         if (rand() < p)
2421                 {
2422                 path = select_cloud_model("Stratus (structured)","large");
2423                 compat_layer.create_cloud(path, lat, lon, balt, 0.0);
2424                 counter = 0;
2425                 i = i+1;
2426                 }
2427         
2428         } # end while
2429
2430 }
2431
2432
2433 ###########################################################
2434 # place a cloud streak 
2435 ###########################################################
2436
2437 var create_streak = func (type, blat, blong, balt, alt_var, nx, xoffset, edgex, x_var, ny, yoffset, edgey, y_var, direction, tri) {
2438
2439 var flag = 0;
2440 var path = "Models/Weather/blank.ac";
2441 calc_geo(blat);
2442 var dir = direction * math.pi/180.0;
2443
2444 var ymin = -0.5 * ny * yoffset;
2445 var xmin = -0.5 * nx * xoffset;
2446 var xinc = xoffset * (tri-1.0) /ny;
2447  
2448 var jlow = int(nx*edgex);
2449 var ilow = int(ny*edgey);
2450
2451
2452 for (var i=0; i<ny; i=i+1)
2453         {
2454         var y = ymin + i * yoffset; 
2455         
2456         for (var j=0; j<nx; j=j+1)
2457                 {
2458                 var y0 = y + y_var * 2.0 * (rand() -0.5);
2459                 var x = xmin + j * (xoffset + i * xinc) + x_var * 2.0 * (rand() -0.5);
2460                 var lat = blat + m_to_lat * (y0 * math.cos(dir) - x * math.sin(dir));
2461                 var long = blong + m_to_lon * (x * math.cos(dir) + y0 * math.sin(dir));
2462
2463                 var alt = balt + alt_var * 2 * (rand() - 0.5);
2464                 
2465                 flag = 0;
2466                 var rn = 6.0 * rand();
2467
2468                 if (((j<jlow) or (j>(nx-jlow-1))) and ((i<ilow) or (i>(ny-ilow-1)))) # select a small or no cloud               
2469                         {
2470                         if (rn > 2.0) {flag = 1;} else {path = select_cloud_model(type,"small");}
2471                         }
2472                 if ((j<jlow) or (j>(nx-jlow-1)) or (i<ilow) or (i>(ny-ilow-1)))         
2473                         {
2474                         if (rn > 5.0) {flag = 1;} else {path = select_cloud_model(type,"small");}
2475                         }
2476                 else    { # select a large cloud
2477                         if (rn > 5.0) {flag = 1;} else {path = select_cloud_model(type,"large");}
2478                         }
2479
2480
2481                 if (flag==0){
2482                         if (thread_flag == 1)
2483                                 {create_cloud_vec(path, lat, long, alt, 0.0);}
2484                         else
2485                                 {compat_layer.create_cloud(path, lat, long, alt, 0.0);}
2486                         
2487
2488                                 }
2489                 }
2490
2491         } 
2492
2493 }
2494
2495
2496
2497
2498
2499
2500
2501 ###########################################################
2502 # place a cloud layer with a gap in the middle
2503 # (useful to reduce cloud count in large thunderstorms)
2504 ###########################################################
2505
2506 var create_hollow_layer = func (type, blat, blon, balt, bthick, rx, ry, phi, density, edge, gap_fraction) {
2507
2508
2509 var i = 0;
2510 var area = math.pi * rx * ry;
2511 var n = int(area/80000000.0 * 100 * density);
2512 var path = "Models/Weather/blank.ac";
2513
2514 phi = phi * math.pi/180.0;
2515
2516 if (contains(cloud_vertical_size_map, type)) 
2517                 {var alt_offset = cloud_vertical_size_map[type]/2.0 * m_to_ft;}
2518         else {var alt_offset = 0.0;}
2519
2520 while(i<n)
2521         {
2522         var x = rx * (2.0 * rand() - 1.0); 
2523         var y = ry * (2.0 * rand() - 1.0); 
2524         var alt = balt + bthick * rand() + 0.8 * alt_offset;
2525         var res = (x*x)/(rx*rx) + (y*y)/(ry*ry);
2526         
2527
2528         if ((res < 1.0) and (res > (gap_fraction * gap_fraction)))
2529                 {
2530                 var lat = blat + m_to_lat * (y * math.cos(phi) - x * math.sin(phi));
2531                 var lon = blon + m_to_lon * (x * math.cos(phi) + y * math.sin(phi));
2532                 if (res > ((1.0 - edge) * (1.0- edge)))
2533                         {
2534                         if (rand() > 0.4) {
2535                                 path = select_cloud_model(type,"small");
2536                                 compat_layer.create_cloud(path, lat, lon, alt, 0.0);
2537                                 }
2538                         }
2539                 else {
2540                         path = select_cloud_model(type,"large");
2541                         if (thread_flag == 1)
2542                                 {create_cloud_vec(path, lat, lon, alt, 0.0);}
2543                         else 
2544                                 {compat_layer.create_cloud(path, lat, lon, alt, 0.0);}
2545                         }
2546                 i = i + 1;
2547                 }
2548         else    # we are in the central gap region
2549                 {
2550                 i = i + 1;
2551                 }
2552         }
2553
2554 i = 0;
2555
2556
2557 }
2558
2559
2560
2561
2562 ###########################################################
2563 # place a cloud box
2564 ###########################################################
2565
2566
2567 var create_cloudbox = func (type, blat, blon, balt, dx,dy,dz,n, f_core, r_core, h_core, n_core, f_bottom, h_bottom, n_bottom) {
2568
2569 var phi = 0;
2570
2571 # first get core coordinates
2572
2573 var core_dx = dx * f_core;
2574 var core_dy = dy * f_core;
2575 var core_dz = dz * h_core;
2576
2577 var core_x_offset = (1.0 * rand() - 0.5) *  ((dx - core_dx) * r_core);
2578 var core_y_offset = (1.0 * rand() - 0.5) *  ((dy - core_dy) * r_core);
2579
2580 # get the bottom geometry
2581
2582 var bottom_dx = dx * f_bottom;
2583 var bottom_dy = dy * f_bottom;
2584 var bottom_dz = dz * h_bottom;
2585
2586 var bottom_offset = 400.0; # in practice, need a small shift
2587
2588 # fill the main body of the box
2589
2590 for (var i=0; i<n; i=i+1)
2591         {
2592
2593         var x = 0.5 * dx * (2.0 * rand() - 1.0); 
2594         var y = 0.5 * dy * (2.0 * rand() - 1.0);
2595         
2596         # veto in core region
2597         if ((x > core_x_offset - 0.5 * core_dx) and (x < core_x_offset + 0.5 * core_dx))
2598                 {
2599                 if ((y > core_y_offset - 0.5 * core_dy) and (y < core_y_offset + 0.5 * core_dy))
2600                         {
2601                         i = i -1;
2602                         continue;
2603                         }
2604                 }
2605          
2606         var alt = balt + bottom_dz + bottom_offset +  dz * rand();
2607         
2608         var lat = blat + m_to_lat * (y * math.cos(phi) - x * math.sin(phi));
2609         var lon = blon + m_to_lon * (x * math.cos(phi) + y * math.sin(phi));
2610
2611         var path = select_cloud_model(type,"standard");
2612
2613         
2614         create_cloud_vec(path, lat, lon, alt, 0.0);
2615
2616         }
2617
2618 # fill the core region
2619
2620 for (var i=0; i<n_core; i=i+1)
2621         {
2622         var x = 0.5 * core_dx * (2.0 * rand() - 1.0); 
2623         var y = 0.5 * core_dy * (2.0 * rand() - 1.0);
2624         var alt = balt + bottom_dz + bottom_offset + core_dz * rand();
2625
2626
2627         var lat = blat + m_to_lat * (y * math.cos(phi) - x * math.sin(phi));
2628         var lon = blon + m_to_lon * (x * math.cos(phi) + y * math.sin(phi));
2629
2630         var path = select_cloud_model(type,"core");
2631
2632         if (thread_flag == 1)
2633                         {create_cloud_vec(path, lat, lon, alt, 0.0);}
2634                 else
2635                         {compat_layer.create_cloud(path, lat, lon, alt, 0.0);}
2636
2637         }
2638
2639 # fill the bottom region
2640
2641
2642 for (var i=0; i<n_bottom; i=i+1)
2643         {
2644         var x = 0.5 * bottom_dx * (2.0 * rand() - 1.0); 
2645         var y = 0.5 * bottom_dy * (2.0 * rand() - 1.0);
2646         var alt = balt + bottom_dz * rand();
2647
2648
2649         var lat = blat + m_to_lat * (y * math.cos(phi) - x * math.sin(phi));
2650         var lon = blon + m_to_lon * (x * math.cos(phi) + y * math.sin(phi));
2651
2652         var path = select_cloud_model(type,"bottom");
2653
2654         if (thread_flag == 1)
2655                         {create_cloud_vec(path, lat, lon, alt, 0.0);}
2656                 else
2657                         {compat_layer.create_cloud(path, lat, lon, alt, 0.0);}
2658
2659         }
2660
2661
2662 }
2663
2664
2665
2666 ###########################################################
2667 # terrain presampling initialization
2668 ###########################################################
2669
2670 var terrain_presampling_start = func (blat, blon, nc, size, alpha) {
2671
2672 # terrain presampling start is always used the first time, and initializes
2673 # the hard-coded routine if that is available since the hard-coded routine cannot
2674 # be yet read out on startup
2675         
2676 # initialize the result vector
2677
2678 setsize(terrain_n,40);
2679 for(var j=0;j<40;j=j+1){terrain_n[j]=0;}
2680
2681 if (thread_flag == 1)
2682         {
2683         var status = getprop(lw~"tmp/presampling-status");
2684         if (status != "idle") # we try a second later
2685                 {
2686                 settimer( func {terrain_presampling_start(blat, blon, nc, size, alpha);},1.00);
2687                 return;
2688                 }
2689         else    
2690                 {
2691                 setprop(lw~"tmp/presampling-status", "sampling");
2692                 terrain_presampling_loop (blat, blon, nc, size, alpha);
2693                 }
2694         }
2695 else
2696         {
2697         terrain_presampling(blat, blon, nc, size, alpha);
2698         terrain_presampling_analysis();
2699         setprop(lw~"tmp/presampling-status", "finished");
2700         }
2701         
2702 if (compat_layer.features.terrain_presampling == 1)
2703         {
2704         print("Starting hard-coded terrain presampling");
2705         setprop("/environment/terrain/area[0]/enabled",1);
2706         setprop(lw~"tmp/presampling-status", "sampling");
2707         setprop("/environment/terrain/area[0]/enabled", 1 );
2708         setprop("/environment/terrain/area[0]/input/latitude-deg", blat );
2709         setprop("/environment/terrain/area[0]/input/longitude-deg", blon );
2710         setprop("/environment/terrain/area[0]/input/use-aircraft-position",1);
2711         setprop("/environment/terrain/area[0]/input/radius-m",45000.0);
2712
2713         setprop("/environment/terrain/area[0]/output/valid", 0 );
2714
2715         }
2716 }
2717
2718 ###########################################################
2719 # terrain presampling loop
2720 ###########################################################
2721
2722 var terrain_presampling_loop = func (blat, blon, nc, size, alpha) {
2723
2724 if ((local_weather_running_flag == 0) and (local_weather_startup_flag == 0)) {return;}
2725
2726
2727 var n = 25;
2728 var n_out = 25;
2729 if (local_weather.features.fast_geodinfo == 0)
2730         {
2731         # dynamically drop accuracy if framerate is low
2732
2733         var dt = getprop("/sim/time/delta-sec");
2734
2735         if (dt > 0.2) # we have below 20 fps
2736                 {n = 5;}
2737         else if (dt > 0.1) # we have below 10 fps
2738                 {n = 10;}
2739         else if (dt > 0.05) # we have below 5 fps
2740                 {n = 15;}
2741         }
2742 else
2743         {
2744         n = 250; n_out = 250;
2745         }
2746
2747 if (nc <= 0) # we're done and may analyze the result
2748         {
2749         terrain_presampling_analysis();
2750         if (debug_output_flag == 1) 
2751                 {print("Presampling done!");}
2752         setprop(lw~"tmp/presampling-status", "finished");
2753         return;
2754         }
2755
2756 terrain_presampling(blat, blon, n, size, alpha);
2757
2758 settimer( func {terrain_presampling_loop(blat, blon, nc-n_out, size, alpha) },0);
2759 }
2760
2761
2762 ###########################################################
2763 # terrain presampling routine
2764 ###########################################################
2765
2766 var terrain_presampling = func (blat, blon, ntries, size, alpha) {
2767
2768 var phi = alpha * math.pi/180.0;
2769 var elevation = 0.0;
2770
2771 var lat_vec = [];
2772 var lon_vec = [];
2773 var lat_lon_vec = [];
2774
2775
2776 for (var i=0; i<ntries; i=i+1)
2777         {
2778         var x = (2.0 * rand() - 1.0) * size;
2779         var y = (2.0 * rand() - 1.0) * size; 
2780         
2781         append(lat_vec, blat + (y * math.cos(phi) - x * math.sin(phi)) * m_to_lat);
2782         append(lon_vec, blon + (x * math.cos(phi) + y * math.sin(phi)) * m_to_lon);
2783         }
2784         
2785         
2786 var elevation_vec = compat_layer.get_elevation_array(lat_vec, lon_vec);
2787         
2788         
2789 for (var i=0; i<ntries;i=i+1)
2790         {
2791         for(var j=0;j<30;j=j+1)
2792                 {
2793                 if ((elevation_vec[i] != -1.0) and (elevation_vec[i] < 500.0 * (j+1))) 
2794                         {terrain_n[j] = terrain_n[j]+1;  break;}
2795                 }
2796                 
2797         }
2798
2799
2800
2801 }
2802
2803 ###########################################################
2804 # terrain presampling analysis
2805 ###########################################################
2806
2807 var terrain_presampling_analysis = func {
2808
2809 if ((compat_layer.features.terrain_presampling_active == 0) or (getprop(lw~"tiles/tile-counter") == 0))
2810         {
2811         var sum = 0;
2812         var alt_mean = 0;
2813         var alt_med = 0;
2814         var alt_20 = 0;
2815         var alt_min = 0;
2816         var alt_low_min = 0;
2817
2818
2819         for (var i=0; i<40;i=i+1)
2820                 {sum = sum + terrain_n[i];}
2821
2822         var n_tot = sum;
2823
2824         sum = 0;
2825         for (var i=0; i<40;i=i+1)
2826                 {
2827                 sum = sum + terrain_n[i];
2828                 if (sum > int(0.5 *n_tot)) {alt_med = i * 500.0; break;}                
2829                 }
2830
2831         sum = 0;
2832         for (var i=0; i<40;i=i+1)
2833                 {
2834                 sum = sum + terrain_n[i];
2835                 if (sum > int(0.3 *n_tot)) {alt_20 = i * 500.0; break;}         
2836                 }
2837
2838
2839         for (var i=0; i<40;i=i+1) {alt_mean = alt_mean + terrain_n[i] * i * 500.0;}
2840         alt_mean = alt_mean/n_tot;
2841
2842         for (var i=0; i<40;i=i+1) {if (terrain_n[i] > 0) {alt_min = i * 500.0; break;}}
2843
2844         var n_max = 0;
2845         sum = 0;
2846
2847         for (var i=0; i<39;i=i+1) 
2848                 {
2849                 sum = sum + terrain_n[i];
2850                 if (terrain_n[i] > n_max) {n_max = terrain_n[i];}
2851                 if ((n_max > terrain_n[i+1]) and (sum > int(0.3*n_tot)))
2852                         {alt_low_min = i * 500; break;}
2853                 }
2854         }
2855 else
2856         {
2857 #       print("Hard-coded sampling...");
2858         var n_tot = getprop("/environment/terrain/area[0]/input/max-samples");
2859         var alt_mean = getprop("/environment/terrain/area[0]/output/alt-mean-ft");
2860         var alt_med = getprop("/environment/terrain/area[0]/output/alt-median-ft");
2861         var alt_min = getprop("/environment/terrain/area[0]/output/alt-min-ft");
2862         var alt_20 = getprop("/environment/terrain/area[0]/output/alt-offset-ft");
2863         }
2864
2865 if (debug_output_flag == 1) 
2866         {print("Terrain presampling analysis results:");
2867         print("total: ",n_tot," mean: ",alt_mean," median: ",alt_med," min: ",alt_min, " alt_20: ", alt_20);}
2868
2869
2870
2871 setprop(lw~"tmp/tile-alt-offset-ft",alt_20);
2872 setprop(lw~"tmp/tile-alt-median-ft",alt_med);
2873 setprop(lw~"tmp/tile-alt-min-ft",alt_min);
2874 setprop(lw~"tmp/tile-alt-mean-ft",alt_mean);
2875 setprop(lw~"tmp/tile-alt-layered-ft",0.5 * (alt_min + alt_20));
2876
2877 append(alt_50_array, alt_med);
2878 append(alt_20_array, alt_20); 
2879 append(alt_min_array, alt_min);
2880 append(alt_mean_array, alt_mean);
2881
2882
2883 current_mean_alt = 0.5 * (current_mean_alt + alt_20);
2884
2885
2886 }
2887
2888
2889
2890 ###########################################################
2891 # wave conditions search
2892 ###########################################################
2893
2894 var wave_detection_loop = func (blat, blon, nx, alpha) {
2895
2896 if (local_weather_running_flag == 0) {return;}
2897
2898 var phi = alpha * math.pi/180.0;
2899 var elevation = 0.0;
2900 var ny = 20;
2901
2902
2903 for (var i=0; i<ny; i=i+1)
2904         {
2905         var x = 5000.0;
2906         var y = -20000.0 + i * 2000.0;
2907         
2908         var lat = blat + (y * math.cos(phi) - x * math.sin(phi)) * m_to_lat;
2909         var lon = blon + (x * math.cos(phi) + y * math.sin(phi)) * m_to_lon;
2910
2911         elevation = compat_layer.get_elevation(lat, lon);
2912
2913         print(elevation);       
2914
2915         }
2916
2917
2918 }
2919
2920 ###########################################################
2921 # detailed altitude determination for convective calls
2922 # clouds follow the terrain to some degree, but not excessively so
2923 ###########################################################
2924
2925 var get_convective_altitude = func (balt, elevation, tile_index, grad) {
2926
2927
2928 var alt_offset = alt_20_array[tile_index - 1];
2929 var alt_median = alt_50_array[tile_index - 1];
2930
2931 # get the maximal shift
2932 var alt_variation = alt_median - alt_offset;
2933
2934 # always get some amount of leeway
2935 if (alt_variation < 500.0) {alt_variation = 500.0;}
2936
2937 # get the correction to the maximal shift by detailed terrain
2938
2939 if (detailed_terrain_interaction_flag == 1)
2940         {
2941         var gradfact = get_gradient_factor(grad);
2942         
2943         if ((local_weather.wind_model_flag == 1) or (local_weather.wind_model_flag == 3))
2944                 {
2945                 var windspeed = tile_wind_speed[0];
2946                 }
2947         else if ((local_weather.wind_model_flag ==2) or (local_weather.wind_model_flag == 4) or (local_weather.wind_model_flag == 5))
2948                 {
2949                 var windspeed = tile_wind_speed[tile_index-1];
2950                 }
2951
2952         var gradfact = ((gradfact - 1.0)  * windspeed) + 1.0;
2953         #print("gradfact: ", gradfact);
2954         }
2955 else
2956         {
2957         var gradfact = 1.0;
2958         }
2959
2960 var alt_variation = alt_variation * gradfact;
2961
2962 # get the difference between offset and foot point
2963 var alt_diff = elevation - alt_offset;
2964
2965 # now get the elevation-induced shift
2966
2967 var fraction = alt_diff / alt_variation;
2968
2969 if (fraction > 1.0) {fraction = 1.0;} # no placement above maximum shift
2970 if (fraction < 0.0) {fraction = 0.0;} # no downward shift
2971
2972 # get the cloud base
2973
2974 var cloudbase = balt - alt_offset;
2975
2976 var alt_above_terrain = balt - elevation;
2977
2978 # the shift strength is weakened if the layer is high above base elevation
2979 # the reference altitude is 1000 ft, anything higher has less sensitivity to terrain
2980
2981 var shift_strength = 1000.0/alt_above_terrain; 
2982
2983 if (shift_strength > 1.0) {shift_strength = 1.0;} # no enhancement for very low layers 
2984 if (shift_strength < 0.0) {shift_strength = 1.0;} # this shouldn't happen, but just in case...
2985
2986 if (alt_diff > alt_variation) {alt_diff = alt_variation;} # maximal shift is given by alt_variation
2987
2988 # print("balt: ", balt, " new alt: ", balt + shift_strength * alt_diff * fraction);
2989
2990 return balt + shift_strength * alt_diff * fraction;
2991
2992 }
2993
2994
2995 ###########################################################
2996 # detailed terrain gradient determination in wind direction
2997 ###########################################################
2998
2999
3000 var get_terrain_gradient = func (lat, lon, elevation1, phi, dist) {
3001
3002
3003 # get the first elevation
3004 # var elevation1 = compat_layer.get_elevation(lat,lon);
3005
3006 # look <dist> upwind to learn about the history of the cloud
3007 var elevation2 = compat_layer.get_elevation(lat+weather_tiles.get_lat(0.0,dist,phi), lon+weather_tiles.get_lon(0.0,dist,phi));
3008
3009 return (elevation2 - elevation1)/(dist * m_to_ft);
3010 }
3011
3012 ###########################################################
3013 # enhancement of the placement altitude due to terrain
3014 ###########################################################
3015
3016 var get_gradient_factor = func (grad) {
3017
3018 if (grad > 0.0)
3019         {return 1.0;}
3020 else
3021         {
3022         return 1.0 -2.0 * grad;
3023         }
3024 }
3025
3026
3027 ###########################################################
3028 # suppression of placement in lee terrain
3029 ###########################################################
3030
3031 var get_lee_bias = func (grad, tile_index) {
3032
3033
3034 if ((local_weather.wind_model_flag == 1) or (local_weather.wind_model_flag == 3))
3035                 {
3036                 var windspeed = tile_wind_speed[0];
3037                 }
3038         else if ((local_weather.wind_model_flag ==2) or (local_weather.wind_model_flag == 4) or (local_weather.wind_model_flag == 5))
3039                 {
3040                 var windspeed = tile_wind_speed[tile_index-1];
3041                 }
3042
3043
3044 if (grad < 0.0)
3045         {return 1.0;}
3046 else
3047         {
3048         var lee_bias = 1.0 - (grad * 0.2 * windspeed);
3049         }
3050 if (lee_bias < 0.2) {lee_bias = 0.2;}
3051
3052 return lee_bias;
3053 }
3054
3055 ###########################################################
3056 # enhancement of Cumulus in above average altitude
3057 ###########################################################
3058
3059
3060 var get_terrain_altitude_factor = func (tile_index, balt, elevation) {
3061
3062
3063 var alt_mean = alt_mean_array[tile_index -1];
3064 var alt_base = alt_20_array[tile_index -1];
3065
3066 var alt_layer = balt - alt_base;
3067 var alt_above_terrain = balt - elevation;
3068 var alt_above_mean = balt - alt_mean;
3069
3070 # the cloud may still be above terrain even if the layer altitude is negative, but we want to avoid neg. factors here
3071
3072 if (alt_above_terrain < 0.0) {alt_above_terrain = 0.0;}
3073
3074 var norm_alt_diff = (alt_above_mean - alt_above_terrain)/alt_layer;
3075
3076 if (norm_alt_diff > 0.0)
3077                 {
3078                 var terrain_altitude_factor = 1.0 + 2.0 * norm_alt_diff;
3079                 }
3080         else
3081                 {
3082                 var terrain_altitude_factor = 1.0/(1.0 - 5.0 * norm_alt_diff);
3083                 }
3084
3085 if (terrain_altitude_factor > 3.0) {terrain_altitude_factor = 3.0;}
3086 if (terrain_altitude_factor < 0.1) {terrain_altitude_factor = 0.1;}
3087
3088 return terrain_altitude_factor;
3089 }
3090
3091
3092 var get_terrain_strength_factor = func (terrain_altitude_factor) {
3093
3094 return  1.0+ (0.5 * (terrain_altitude_factor-1.0));
3095
3096 }
3097
3098
3099 ###########################################################
3100 # terrain presampling listener dispatcher
3101 ###########################################################
3102
3103 var manage_presampling = func {
3104
3105
3106
3107 var status = getprop(lw~"tmp/presampling-status");
3108
3109
3110 # we only take action when the analysis is done
3111 if (status != "finished") {return;} 
3112
3113 if (getprop(lw~"tiles/tile-counter") == 0) # we deal with a tile setup call from the menu
3114         {
3115         set_tile();
3116         }
3117 else    # the tile setup call came from weather_tile_management
3118         {
3119         var lat = getprop(lw~"tiles/tmp/latitude-deg");
3120         var lon = getprop(lw~"tiles/tmp/longitude-deg");
3121         var code = getprop(lw~"tiles/tmp/code");
3122         var dir_index = getprop(lw~"tiles/tmp/dir-index");      
3123
3124         weather_tile_management.generate_tile(code, lat, lon, dir_index);
3125         }
3126
3127
3128 # set status to idle again
3129
3130 setprop(lw~"tmp/presampling-status", "idle");
3131
3132 }
3133
3134
3135 ###########################################################
3136 # hardcoded terrain presampling listener dispatcher
3137 ###########################################################
3138
3139 var manage_hardcoded_presampling = func {
3140
3141 var status = getprop("/environment/terrain/area[0]/enabled");
3142
3143 print("Hard-coded terrain presampling status: ", status);
3144
3145 # no action unless the sampler has finished
3146 if (status ==0) {return;}
3147
3148 # no action if the sampler hasn't been started
3149
3150 if (getprop(lw~"tmp/presampling-status") != "sampling") {return;}
3151
3152 terrain_presampling_analysis();
3153 if (debug_output_flag == 1) 
3154                 {print("Presampling done!");}
3155 setprop(lw~"tmp/presampling-status", "finished");
3156
3157
3158 }
3159
3160 ###########################################################
3161 # set wind model flag
3162 ###########################################################
3163
3164 var set_wind_model_flag = func {
3165
3166 var wind_model = getprop(lw~"config/wind-model");
3167
3168 if (wind_model == "constant") {wind_model_flag = 1;}
3169 else if (wind_model == "constant in tile") {wind_model_flag =2;}
3170 else if (wind_model == "aloft interpolated") {wind_model_flag =3; }
3171 else if (wind_model == "airmass interpolated") {wind_model_flag =4;}
3172 else if (wind_model == "aloft waypoints") {wind_model_flag =5;}
3173 else {print("Wind model not implemented!"); wind_model_flag =1;}
3174
3175
3176 }
3177
3178
3179 ###########################################################
3180 # set texture mix for convective clouds
3181 ###########################################################
3182
3183 var set_texture_mix = func {
3184
3185 var thermal_properties = getprop(lw~"config/thermal-properties");
3186 thermal_conditions = thermal_properties;
3187
3188 convective_texture_mix = -(thermal_properties - 1.0) * 0.4;
3189
3190 if (convective_texture_mix < -0.2) {convective_texture_mix = -0.2;}
3191 if (convective_texture_mix > 0.2) {convective_texture_mix = 0.2;}
3192
3193 lowest_layer_turbulence = 0.7 - thermal_properties;
3194 if (lowest_layer_turbulence < 0.0) {lowest_layer_turbulence = 0.0;}
3195 }
3196
3197 ###########################################################
3198 # create an effect volume
3199 ###########################################################
3200
3201 var create_effect_volume = func (geometry, lat, lon, r1, r2, phi, alt_low, alt_high, vis, rain, snow, turb, lift, lift_flag, sat) {
3202
3203
3204 var ev = effectVolume.new (geometry, lat, lon, r1, r2, phi, alt_low, alt_high, vis, rain, snow, turb, lift, lift_flag, sat);
3205 ev.index = getprop(lw~"tiles/tile-counter");
3206 ev.active_flag = 0;
3207
3208
3209 if (vis < 0.0) {ev.vis_flag = 0;} else {ev.vis_flag = 1;}
3210 if (rain < 0.0) {ev.rain_flag = 0;} else {ev.rain_flag = 1;}
3211 if (snow < 0.0) {ev.snow_flag = 0;} else {ev.snow_flag = 1;}
3212 if (turb < 0.0) {ev.turb_flag = 0;} else {ev.turb_flag = 1;}
3213 if (lift_flag ==  0.0) {ev.lift_flag = 0;} else {ev.lift_flag = 1;}
3214 if (sat < 0.0) {ev.sat_flag = 0;} else {ev.sat_flag = 1;}
3215 if (sat > 1.0) {sat = 1.0;}
3216
3217 if (lift_flag == -2) # we create a thermal by function
3218         {
3219         ev.lift_flag = 2;
3220         ev.radius = 0.8 * r1;
3221         ev.height = alt_high * 0.87;
3222         ev.cn = 0.7 + rand() * 0.2;
3223         ev.sh = 0.7 + rand() * 0.2;
3224         ev.max_lift = lift;
3225         ev.f_lift_radius = 0.7 + rand() * 0.2;
3226         if (dynamics_flag == 1) # globals set by the convective system
3227                 {
3228                 ev.flt = cloud_fractional_lifetime;
3229                 ev.evolution_timestamp = cloud_evolution_timestamp;
3230                 }
3231         }
3232
3233 if (lift_flag == -3) # we create a wave lift
3234         {
3235         ev.lift_flag = 3;
3236         ev.height = 10000.0; # scale height in ft
3237         ev.max_lift = lift;
3238         ev.index = 0; # static objects are assigned tile id zero
3239         }
3240
3241 # set a timestamp if needed
3242
3243 if (dynamics_flag == 1)
3244         {
3245         ev.timestamp = weather_dynamics.time_lw;
3246         }
3247
3248 # and add to the counter
3249 setprop(lw~"effect-volumes/number",getprop(lw~"effect-volumes/number")+1);
3250
3251 append(effectVolumeArray,ev);
3252 }
3253
3254
3255
3256
3257
3258 ###########################################################
3259 # set a weather station for interpolation
3260 ###########################################################
3261
3262 var set_weather_station = func (lat, lon, alt, vis, T, D, p) {
3263
3264 var s = weatherStation.new (lat, lon, alt, vis, T, D, p);
3265 s.index = getprop(lw~"tiles/tile-counter");
3266 s.weight = 0.02;
3267
3268 # set a timestamp if needed
3269