Attachment 'hc08t.c'
Download 1 // hc08t.c
2 // compute heating and cooling, with timestamp and color scale
3 //
4 // gcc -o hc08t hc08t.c -lgd -lpng -lm
5
6 #define NAME "hc08t"
7
8 #define RATE 12 // Flash frame rate ( 1560x )
9 #define TIME 1800 // Total time
10 #define NPLOT 360 // Number of plots
11 #define NCOMP 70 // Computation passes per plot
12 #define PPLOT 2400 // Plot width in KM (may be ugly!)
13 #define OPLOT -200 // Left side of plot
14 #define TOFF 29 // First payload start time
15
16 #define TEMP0 340.0 // K Starting temperature
17 #define TEMPA 250.0 // K Ambient temperature
18 #define PV 10800.0 // m/s Payload maximum velocity
19 #define PA 30.0 // m/s2 Payload acceleration
20 #define PT 360.0 // s Payload acceleration time
21 #define PM 5000.0 // kg Payload mass
22 #define PL 1944.0 // km Payload acceleration distance
23 #define PS 45.0 // s Payload spacing in time
24 #define PN 30 // Number of computed payloads
25
26 #define RL 5600 // km Rotor Length
27 #define RM 3.0 // kg/m Rotor mass/length
28 #define RR 0.025 // m Rotor radius in meters
29 #define RV 14000.0 // m/s Rotor velocity
30
31 #define EM 0.8 // Emissivity
32 #define SB 5.67E-8 // W/m2K4 Blackbody radiation
33 #define HC 600 // J/kg-K Rotor heat capacity
34 #define PI 3.14159265
35
36 #define MAXTEMP (300+3*255)
37 #define FNT "DejaVuMonoSans"
38 #define FSZ 20
39
40 // MAJOR KLUDGE WARNING! These plot corners depend on what Gnuplot draws:
41 #define XGP_UL 168
42 #define YGP_UL 83
43 #define XGP_LR 974
44 #define YGP_LR 683
45 #define TGP_LO 200 // low temperature
46 #define TGP_HI 1000 // high temperature, must be MAXTEMP
47 #define YHBAND 3
48 #define YHSTEP 5
49
50 #include <gd.h>
51 #include <math.h>
52 #include <string.h>
53 #include <stdio.h>
54 #include <stdlib.h>
55
56 double rtemp[RL] ; // rotor temperature K
57
58 //-----------------------------------------------------------------
59 // Rotor loop modulus function
60
61 int rmax = RL * ( 500000000 / RL ) ; // large offset multiple of RL
62
63 int rmod( int arg ) {
64 return ( ( rmax + arg ) % RL ) ;
65 }
66
67 //-----------------------------------------------------------------
68 int x_to_i ( double xarg ) {
69 return (int) ( XGP_UL + ( (double) ( XGP_LR - YGP_UL ) )
70 * ( (double) ( xarg - OPLOT ) )
71 / ( (double) ( PPLOT - OPLOT ) )
72 ) ;
73 }
74
75 //-----------------------------------------------------------------
76 int main() {
77
78 FILE *datfile ; //
79 char rmmkdir[200] ; //
80 char png2swf[200] ; //
81 char gnuplot[200] ; //
82 char filename[80] ; //
83 char plotlabel[80] ; // plot label
84 char colortempr[MAXTEMP] ; //
85
86 double ptime ; // acceleration time of payload
87 double vp ; // velocity of payload
88 double xp ; // position of payload
89
90 int ipoint ; // point index for plotting, cooling
91 int jpoint ; // point index for plotting, cooling
92 int index = 0 ; // movement of rotor since t0
93 int iplot ; // plot number
94 int icomp ; // computation index
95 int ipay ; // payload index
96 double rtp ; // temporary rtemp
97 double time ; // simulating time
98 double tdelta = 1000.0/RV ; // timestep in seconds
99 double xpoint ; //
100
101 double xpay[PN] ; // viewable payload position
102 double xht0[PN] ; // front of heat wave
103 double xht1[PN] ; // index of payload
104 int cht[PN] ; // heat color
105 double hpass ; // heat pass
106 int npay ; // number of viewable payloads
107 double rl = (double) RL ; // loop roundtrip length
108
109 double dtp = tdelta // payload heating factor
110 * PA * PM // payload force
111 * 0.001 // kilometer of mass
112 / ( RM * HC ) ; // rotor heat capacity J/K
113
114 double heat ; // added heat
115 int ip ; // integer position of payload
116 int ip0 ; // payload position 0, rotor-coords
117 int ip1 ; // payload position 1, rotor-coords
118 double xpi ; // integer part of payload position
119 double xpf ; // fractional part of position
120
121 double dt0 = tdelta // rotor cooling factor
122 * EM * SB // blackbody radiation
123 * 2.0 * PI * RR // surface area
124 / ( RM * HC ) ; // rotor heat capacity J/K
125
126 double dt1 = dt0 // rotor ambient
127 * TEMPA*TEMPA //
128 * TEMPA*TEMPA ; //
129
130 int speedup = (RATE*TIME)/NPLOT ;
131
132 int i, j ; // general index
133 int ix, iy ; // general index for color
134 int hx ; // end of heat wave
135 int itempr ; // temporary integer temperature
136
137 double tangle ; // turnaround angle
138 double tsmooth ; // turnaround angle smoothing
139 int x1,y1,x2,y2 ; // drawing points for turnaround
140
141 gdPoint payload[3] ; // payload triangle
142 gdImagePtr im ; // working image
143
144 FILE *pngin; ; // output file
145 FILE *pngout; ; // output file
146 int black, white, blue, red ; // colors
147
148 int dxpath = (0.5*PA*PT-RV)*PT; // after acceleration
149
150 // --------------------------------------------------------------------
151
152 printf ( "tdelta = %14.4e\n", tdelta );
153 printf ( "dt0 = %14.4e\n", dt0 );
154 printf ( "dt1 = %14.4e\n", dt1 );
155 printf ( "dtp = %14.4e\n", dtp );
156
157 sprintf( rmmkdir, "rm -rf %sdir ; mkdir %sdir", NAME, NAME );
158 system( rmmkdir ) ; // make directory
159
160 // initialize
161 for( ipoint=0 ; ipoint < RL ; ipoint++ ) { rtemp[ipoint] = TEMP0 ; }
162
163 // plot loop
164 for( iplot=0 ; iplot < NPLOT ; iplot++ ) {
165
166 printf( "plot %4d\r", iplot ); fflush( stdout );
167
168 // compute loop per plot
169 for( icomp=0 ; icomp < NCOMP ; icomp++ ) {
170
171 index = iplot*NCOMP + icomp ;
172 time = tdelta * ((double) index ) - TOFF ;
173
174 // cooling
175 for( ipoint=0 ; ipoint < RL ; ipoint++ ) {
176 rtp = rtemp[ ipoint ] ;
177 rtp *= rtp ;
178 rtemp[ ipoint ] -= dt0*rtp*rtp - dt1 ;
179 }
180
181 npay = 0 ;
182 // heating by payloads - start a little later
183 for( ipay = 0 ; ipay < PN ; ipay++ ) {
184 ptime = time - PS * ((double)ipay) ;
185
186 if( ptime < 0.0 ) {
187 xht0[ ipay ] = -1000.0 ;
188 xht1[ ipay ] = -1000.0 ;
189 }
190 else if ( ptime <= PT ) {
191 vp = PA * ptime ;
192 xp = 0.0005 * vp * ptime ;
193 heat = dtp * ( RV - vp ) ;
194 xpi = floor( xp ) ;
195 ip = (int) xpi ;
196 xpf = xp-xpi ;
197 ip0 = rmod( ip - index ) ;
198 ip1 = rmod( ip+1 - index ) ;
199 rtemp[ ip0 ] += heat*(1.0-xpf) ;
200 rtemp[ ip1 ] += heat*( xpf) ;
201 xht0[ ipay ] = xp ;
202 xht1[ ipay ] = 0.001*ptime*RV ; // km
203 cht[ ipay ] = 600 ;
204 xpay[npay++] = xp ;
205 }
206 else { // ptime > PT
207 xht0[ ipay ] = 0.001*(ptime*RV+dxpath) ; // km
208 xht1[ ipay ] = 0.001*ptime*RV ; // km
209 hpass = floor( ( xht1[ ipay ] -OPLOT ) / rl ) ;
210 cht[ ipay ] = 700-100*((int)hpass ) ;
211 if( cht[ ipay ] < 0 ) { cht[ ipay ] = 0 ; }
212 xht0[ ipay ] -= rl*hpass ;
213 xht1[ ipay ] -= rl*hpass ;
214 if ( xht0[ ipay ] < OPLOT ) { xht0[ ipay ] = OPLOT ; }
215 }
216 }
217 } // end of icomp compute loop
218
219 index = (iplot+1)*NCOMP ;
220 datfile = fopen( "tmp.dat", "wb" );
221
222 for( ipoint = 0 ; ipoint < PPLOT ; ipoint++ ) {
223 xpoint = (double)( OPLOT+ipoint) ; // starts before 0
224 ip = rmod( OPLOT+ipoint-index );
225 fprintf( datfile, "%16.4f%16.4f\n", xpoint, rtemp[ ip ] ) ;
226 }
227 fclose( datfile );
228
229 sprintf( gnuplot, "/usr/local/bin/gp %s.cmd", NAME );
230 system( gnuplot );
231
232 // ===================================================================
233 // re-open plot with libgd, add time stamp
234 // This scribbles on top of existing plot, so it is dependent
235 // on UL/LR placement in plot window
236
237 pngin = fopen( "tmp.png", "rb" );
238 im = gdImageCreateFromPng(pngin);
239 fclose( pngin );
240 black = gdImageColorAllocate (im, 0, 0, 0);
241 blue = gdImageColorAllocate (im, 0, 0, 255);
242 red = gdImageColorAllocate (im, 224, 0, 0);
243
244 // color scale ------------------------------------------------------
245 // 0 to 300+3*255 = 1065
246 for( i=0 ; i<300; i++ ) { colortempr[i] = black ; }
247
248 // black to red to yellow to white
249 // only 256 colors for image, assign 153 of them
250 for( i=0 ; i<255; i+=5 ) {
251 colortempr[i+300] = gdImageColorAllocate( im, i, 0, 0 );
252 colortempr[i+555] = gdImageColorAllocate( im, 255, i, 0 );
253 colortempr[i+810] = gdImageColorAllocate( im, 255, 255, i );
254 for( j=1 ; j<5 ; j++ ) {
255 colortempr[i+j+300]=colortempr[i+300];
256 colortempr[i+j+555]=colortempr[i+555];
257 colortempr[i+j+810]=colortempr[i+810];
258 }
259 }
260
261 // right side color temperature scale
262 for( iy=YGP_UL ; iy <= YGP_LR ; iy++ ) {
263 itempr = (int) ( TGP_LO + ( (double) ( TGP_HI - TGP_LO ) )
264 * ( (double) ( iy - YGP_LR ) )
265 / ( (double) ( YGP_UL - YGP_LR ) )
266 ) ;
267
268 if ( itempr < 930 ) {
269 gdImageLine(im, XGP_LR+ 1, iy,
270 XGP_LR+10, iy, colortempr [ itempr ] );
271 }
272 }
273
274 // track temperature
275 for( ix = XGP_UL ; ix <= XGP_LR ; ix++ ) {
276 ipoint = (int) ( ( (double) ( PPLOT - OPLOT ) )
277 * ( (double) ( ix - XGP_UL ) )
278 / ( (double) ( XGP_LR - YGP_UL ) )
279 ) ;
280
281 // forward track
282 // points between (-)OPLOT and PPLOT , minus the index
283
284 ip = rmod( OPLOT+ipoint-index );
285 itempr = (int) ( rtemp[ ip ] );
286 if ( itempr >= MAXTEMP ) { itempr = MAXTEMP-1 ; }
287 gdImageLine(im, ix, YGP_UL+ 7,
288 ix, YGP_UL-10, colortempr[ itempr ] );
289
290 // reverse track
291 // points between RL/2+PPLOT and RL/2-OPLOT, minus the index
292
293 jpoint = RL/2+PPLOT-ipoint ;
294 ip = rmod( OPLOT+jpoint-index );
295 itempr = (int) ( rtemp[ ip ] );
296 if ( itempr >= MAXTEMP ) { itempr = MAXTEMP-1 ; }
297 gdImageLine(im, ix, YGP_UL+43,
298 ix, YGP_UL+50, colortempr[ itempr ] );
299 }
300
301 // turnarounds
302 for( ipoint = PPLOT ; ipoint < ((RL/2)); ipoint++ ) {
303
304 for( tsmooth=0.9 ; tsmooth > -0.01 ; tsmooth -= 0.1 ) { // smoothing
305 tangle = PI* ( tsmooth + (double)(ipoint-PPLOT))
306 / ( (double)((RL/2)-PPLOT)) ;
307
308 // east turnaround
309 ip = rmod( ipoint+OPLOT-index );
310 x1 = XGP_LR +18*sin(tangle);
311 y1 = YGP_UL+25 -18*cos(tangle);
312 x2 = XGP_LR +30*sin(tangle);
313 y2 = YGP_UL+20 -30*cos(tangle);
314 itempr = (int) ( rtemp[ ip ] );
315 gdImageLine(im, x1, y1, x2, y2, colortempr[ itempr ] );
316
317 // west turnaround
318 ip = rmod( RL/2+ipoint+OPLOT-index );
319 x1 = XGP_UL -18*sin(tangle);
320 y1 = YGP_UL+25 +18*cos(tangle);
321 x2 = XGP_UL -30*sin(tangle);
322 y2 = YGP_UL+20 +30*cos(tangle);
323 itempr = (int) ( rtemp[ ip ] );
324 gdImageLine(im, x1, y1, x2, y2, colortempr[ itempr ] );
325 }
326 }
327
328 // payloads on top of forward track
329 for( i = 0 ; i < npay ; i++ ) {
330 ix = x_to_i( xpay[i] );
331
332 payload[0].x = ix ; payload[0].y = YGP_UL - 10 ;
333 payload[1].x = ix-15 ; payload[1].y = YGP_UL - 25 ;
334 payload[2].x = ix+15 ; payload[2].y = YGP_UL - 25 ;
335 gdImageFilledPolygon( im, payload, 3, blue );
336 }
337
338 // payload heat tracks
339 for( i = 0 ; i < PN ; i++ ) {
340 ix = x_to_i( xht0[i] );
341 hx = x_to_i( xht1[i] );
342 if ( hx > 0 ) {
343 if ( ix < XGP_LR ) {
344 if ( hx > XGP_LR ) { hx = XGP_LR ; }
345 iy = YGP_LR - (i+1)*YHSTEP ;
346 gdImageFilledRectangle( im, ix, iy, hx, iy+YHBAND ,
347 colortempr[ cht[i] ] );
348 }
349 }
350 }
351
352 // ------------------------------------------------------------------
353
354 sprintf( plotlabel, "%4d/%4d sec,%3d vehicles,%4dx animation speedup\n",
355 ((int)time), TIME, PN, speedup );
356
357 gdImageStringFT( im, NULL, // imagespace, bounding box
358 black , FNT, FSZ, 0.0, // color, font, fontsize, angle
359 XGP_UL+2, YGP_UL+36,
360 plotlabel ); // x, y, text
361
362 pngout = fopen( "tmpl.png", "wb");
363 gdImagePngEx( im, pngout, 1 );
364 gdImageDestroy(im);
365 fclose(pngout);
366
367 // move temp file to directory
368 sprintf( filename, "%sdir/b%04d.png", NAME, iplot );
369 rename( "tmpl.png", filename );
370 }
371
372 // make flash video
373 sprintf( png2swf, "png2swf -o %s.swf -r%3d %sdir/*.png", NAME, RATE, NAME );
374 system( png2swf );
375
376 return(0);
377 }
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.