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