Attachment 'hc04.c'
Download 1 // hc04.c
2 // compute heating and cooling
3 //
4 // gcc -o hc04 hc04.c -lm
5
6 #define NAME "hc04"
7 #define CMD "hc"
8
9 #define RATE 12 // Flash frame rate ( 60x )
10 #define TIME 3600 // Total time
11 #define NPLOT 720 // Number of plots
12 #define NCOMP 70 // Computation passes per plot
13 #define PPLOT 2400 // Plot width (may be ugly!)
14 #define OPLOT -200 //
15
16 #define TEMP0 315.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
24 #define PS 24.0 // s Payload spacing in time
25 #define PN 15 // Number of computed payloads
26
27 #define RL 5180 // km Rotor Length
28 #define RM 3.0 // kg/m Rotor mass/length
29 #define RR 0.025 // m Rotor radius in meters
30 #define RV 14000.0 // m/s Rotor velocity
31
32 #define EM 0.8 // Emissivity
33 #define SB 5.67E-8 // W/m2K4 Blackbody radiation
34 #define HC 600 // J/kg-K Rotor heat capacity
35 #define PI 3.14159265
36
37 #include <math.h>
38 #include <string.h>
39 #include <stdio.h>
40 #include <stdlib.h>
41
42 double rtemp[RL] ; // rotor temperature K
43
44 // modulus function
45 int rmax = RL * ( 500000000 / RL ) ; // large offset multiple of RL
46 int rmod( int arg ) { return ( ( rmax + arg ) % RL ) ; }
47
48 //-----------------------------------------------------------------
49 int main() {
50
51 FILE *datfile ; //
52 char rmmkdir[200] ; //
53 char png2swf[200] ; //
54 char gnuplot[200] ; //
55 char filename[80] ; //
56
57 double ptime ; // acceleration time of payload
58 double vp ; // velocity of payload
59 double xp ; // position of payload
60
61 int ipoint ; // point index for plotting, cooling
62 int index = 0 ; // movement of rotor since t0
63 int iplot ; // plot number
64 int icomp ; // computation index
65 int ipay ; // payload index
66 double rtp ; // temporary rtemp
67 double time ; // simulating time
68 double tdelta = 1000.0/RV ; // timestep in seconds
69 double xpoint ; //
70
71 double dtp = tdelta // payload heating factor
72 * PA * PM // payload force
73 * 0.001 // kilometer of mass
74 / ( RM * HC ) ; // rotor heat capacity J/K
75
76 double heat ; // added heat
77 int ip ; // integer position of payload
78 int ip0 ; // payload position 0, rotor-coords
79 int ip1 ; // payload position 1, rotor-coords
80 double xpi ; // integer part of payload position
81 double xpf ; // fractional part of position
82
83 double dt0 = tdelta // rotor cooling factor
84 * EM * SB // blackbody radiation
85 * 2.0 * PI * RR // surface area
86 / ( RM * HC ) ; // rotor heat capacity J/K
87
88 double dt1 = dt0 // rotor ambient
89 * TEMPA*TEMPA //
90 * TEMPA*TEMPA ; //
91
92 // --------------------------------------------------------------------
93
94 printf ( "tdelta = %14.4e\n", tdelta );
95 printf ( "dt0 = %14.4e\n", dt0 );
96 printf ( "dt1 = %14.4e\n", dt1 );
97 printf ( "dtp = %14.4e\n", dtp );
98
99 sprintf( rmmkdir, "rm -rf %sdir ; mkdir %sdir", NAME, NAME );
100 system( rmmkdir ) ; // make directory
101
102 // initialize
103 for( ipoint=0 ; ipoint < RL ; ipoint++ ) { rtemp[ipoint] = TEMP0 ; }
104
105 // plot loop
106 for( iplot=0 ; iplot < NPLOT ; iplot++ ) {
107
108 printf( "plot %4d\r", iplot ); fflush( stdout );
109
110 // compute loop per plot
111 for( icomp=0 ; icomp < NCOMP ; icomp++ ) {
112
113 index = iplot*NCOMP + icomp ;
114 time = tdelta * ((double) index );
115
116 // cooling
117 for( ipoint=0 ; ipoint < RL ; ipoint++ ) {
118 rtp = rtemp[ ipoint ] ;
119 rtp *= rtp ;
120 rtemp[ ipoint ] -= dt0*rtp*rtp - dt1 ;
121 }
122
123 // heating by payloads - start a little later
124 for( ipay = 1 ; ipay <= PN ; ipay++ ) {
125 ptime = time - PS * ((double)ipay) ;
126 if( ( ptime >= 0.0 ) && ( ptime <= PT ) ) {
127
128 vp = PA * ptime ;
129 heat = dtp * ( RV - vp ) ;
130
131 xp = 0.0005 * vp * ptime ;
132 xpi = floor( xp ) ;
133 ip = (int) xpi ;
134 xpf = xp-xpi ;
135 ip0 = rmod( ip - index ) ;
136 ip1 = rmod( ip+1 - index ) ;
137
138 rtemp[ ip0 ] += heat*(1.0-xpf) ;
139 rtemp[ ip1 ] += heat*( xpf) ;
140 }
141 }
142 } // end of icomp compute loop
143
144 index = (iplot+1)*NCOMP ;
145 datfile = fopen( "tmp.dat", "wb" );
146
147 for( ipoint = 0 ; ipoint < PPLOT ; ipoint++ ) {
148 xpoint = (double)( OPLOT+ipoint) ; // starts before 0
149 ip = rmod( OPLOT+ipoint-index );
150 fprintf( datfile, "%16.4f%16.4f\n", xpoint, rtemp[ ip ] ) ;
151 }
152 fclose( datfile );
153
154 sprintf( gnuplot, "/usr/local/bin/gp %s.cmd", CMD );
155 system( gnuplot );
156
157 sprintf( filename, "%sdir/b%04d.png", NAME, iplot );
158 rename( "tmp.png", filename );
159 }
160
161 // make flash video
162 sprintf( png2swf, "png2swf -o %s.swf -r%3d %sdir/*.png", NAME, RATE, NAME );
163 system( png2swf );
164
165 return(0);
166 }
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.