// hc08t.c // compute heating and cooling, with timestamp and color scale // // gcc -o hc08t hc08t.c -lgd -lpng -lm #define NAME "hc08t" #define RATE 12 // Flash frame rate ( 1560x ) #define TIME 1800 // Total time #define NPLOT 360 // Number of plots #define NCOMP 70 // Computation passes per plot #define PPLOT 2400 // Plot width in KM (may be ugly!) #define OPLOT -200 // Left side of plot #define TOFF 29 // First payload start time #define TEMP0 340.0 // K Starting temperature #define TEMPA 250.0 // K Ambient temperature #define PV 10800.0 // m/s Payload maximum velocity #define PA 30.0 // m/s2 Payload acceleration #define PT 360.0 // s Payload acceleration time #define PM 5000.0 // kg Payload mass #define PL 1944.0 // km Payload acceleration distance #define PS 45.0 // s Payload spacing in time #define PN 30 // Number of computed payloads #define RL 5600 // km Rotor Length #define RM 3.0 // kg/m Rotor mass/length #define RR 0.025 // m Rotor radius in meters #define RV 14000.0 // m/s Rotor velocity #define EM 0.8 // Emissivity #define SB 5.67E-8 // W/m2K4 Blackbody radiation #define HC 600 // J/kg-K Rotor heat capacity #define PI 3.14159265 #define MAXTEMP (300+3*255) #define FNT "DejaVuMonoSans" #define FSZ 20 // MAJOR KLUDGE WARNING! These plot corners depend on what Gnuplot draws: #define XGP_UL 168 #define YGP_UL 83 #define XGP_LR 974 #define YGP_LR 683 #define TGP_LO 200 // low temperature #define TGP_HI 1000 // high temperature, must be MAXTEMP #define YHBAND 3 #define YHSTEP 5 #include #include #include #include #include double rtemp[RL] ; // rotor temperature K //----------------------------------------------------------------- // Rotor loop modulus function int rmax = RL * ( 500000000 / RL ) ; // large offset multiple of RL int rmod( int arg ) { return ( ( rmax + arg ) % RL ) ; } //----------------------------------------------------------------- int x_to_i ( double xarg ) { return (int) ( XGP_UL + ( (double) ( XGP_LR - YGP_UL ) ) * ( (double) ( xarg - OPLOT ) ) / ( (double) ( PPLOT - OPLOT ) ) ) ; } //----------------------------------------------------------------- int main() { FILE *datfile ; // char rmmkdir[200] ; // char png2swf[200] ; // char gnuplot[200] ; // char filename[80] ; // char plotlabel[80] ; // plot label char colortempr[MAXTEMP] ; // double ptime ; // acceleration time of payload double vp ; // velocity of payload double xp ; // position of payload int ipoint ; // point index for plotting, cooling int jpoint ; // point index for plotting, cooling int index = 0 ; // movement of rotor since t0 int iplot ; // plot number int icomp ; // computation index int ipay ; // payload index double rtp ; // temporary rtemp double time ; // simulating time double tdelta = 1000.0/RV ; // timestep in seconds double xpoint ; // double xpay[PN] ; // viewable payload position double xht0[PN] ; // front of heat wave double xht1[PN] ; // index of payload int cht[PN] ; // heat color double hpass ; // heat pass int npay ; // number of viewable payloads double rl = (double) RL ; // loop roundtrip length double dtp = tdelta // payload heating factor * PA * PM // payload force * 0.001 // kilometer of mass / ( RM * HC ) ; // rotor heat capacity J/K double heat ; // added heat int ip ; // integer position of payload int ip0 ; // payload position 0, rotor-coords int ip1 ; // payload position 1, rotor-coords double xpi ; // integer part of payload position double xpf ; // fractional part of position double dt0 = tdelta // rotor cooling factor * EM * SB // blackbody radiation * 2.0 * PI * RR // surface area / ( RM * HC ) ; // rotor heat capacity J/K double dt1 = dt0 // rotor ambient * TEMPA*TEMPA // * TEMPA*TEMPA ; // int speedup = (RATE*TIME)/NPLOT ; int i, j ; // general index int ix, iy ; // general index for color int hx ; // end of heat wave int itempr ; // temporary integer temperature double tangle ; // turnaround angle double tsmooth ; // turnaround angle smoothing int x1,y1,x2,y2 ; // drawing points for turnaround gdPoint payload[3] ; // payload triangle gdImagePtr im ; // working image FILE *pngin; ; // output file FILE *pngout; ; // output file int black, white, blue, red ; // colors int dxpath = (0.5*PA*PT-RV)*PT; // after acceleration // -------------------------------------------------------------------- printf ( "tdelta = %14.4e\n", tdelta ); printf ( "dt0 = %14.4e\n", dt0 ); printf ( "dt1 = %14.4e\n", dt1 ); printf ( "dtp = %14.4e\n", dtp ); sprintf( rmmkdir, "rm -rf %sdir ; mkdir %sdir", NAME, NAME ); system( rmmkdir ) ; // make directory // initialize for( ipoint=0 ; ipoint < RL ; ipoint++ ) { rtemp[ipoint] = TEMP0 ; } // plot loop for( iplot=0 ; iplot < NPLOT ; iplot++ ) { printf( "plot %4d\r", iplot ); fflush( stdout ); // compute loop per plot for( icomp=0 ; icomp < NCOMP ; icomp++ ) { index = iplot*NCOMP + icomp ; time = tdelta * ((double) index ) - TOFF ; // cooling for( ipoint=0 ; ipoint < RL ; ipoint++ ) { rtp = rtemp[ ipoint ] ; rtp *= rtp ; rtemp[ ipoint ] -= dt0*rtp*rtp - dt1 ; } npay = 0 ; // heating by payloads - start a little later for( ipay = 0 ; ipay < PN ; ipay++ ) { ptime = time - PS * ((double)ipay) ; if( ptime < 0.0 ) { xht0[ ipay ] = -1000.0 ; xht1[ ipay ] = -1000.0 ; } else if ( ptime <= PT ) { vp = PA * ptime ; xp = 0.0005 * vp * ptime ; heat = dtp * ( RV - vp ) ; xpi = floor( xp ) ; ip = (int) xpi ; xpf = xp-xpi ; ip0 = rmod( ip - index ) ; ip1 = rmod( ip+1 - index ) ; rtemp[ ip0 ] += heat*(1.0-xpf) ; rtemp[ ip1 ] += heat*( xpf) ; xht0[ ipay ] = xp ; xht1[ ipay ] = 0.001*ptime*RV ; // km cht[ ipay ] = 600 ; xpay[npay++] = xp ; } else { // ptime > PT xht0[ ipay ] = 0.001*(ptime*RV+dxpath) ; // km xht1[ ipay ] = 0.001*ptime*RV ; // km hpass = floor( ( xht1[ ipay ] -OPLOT ) / rl ) ; cht[ ipay ] = 700-100*((int)hpass ) ; if( cht[ ipay ] < 0 ) { cht[ ipay ] = 0 ; } xht0[ ipay ] -= rl*hpass ; xht1[ ipay ] -= rl*hpass ; if ( xht0[ ipay ] < OPLOT ) { xht0[ ipay ] = OPLOT ; } } } } // end of icomp compute loop index = (iplot+1)*NCOMP ; datfile = fopen( "tmp.dat", "wb" ); for( ipoint = 0 ; ipoint < PPLOT ; ipoint++ ) { xpoint = (double)( OPLOT+ipoint) ; // starts before 0 ip = rmod( OPLOT+ipoint-index ); fprintf( datfile, "%16.4f%16.4f\n", xpoint, rtemp[ ip ] ) ; } fclose( datfile ); sprintf( gnuplot, "/usr/local/bin/gp %s.cmd", NAME ); system( gnuplot ); // =================================================================== // re-open plot with libgd, add time stamp // This scribbles on top of existing plot, so it is dependent // on UL/LR placement in plot window pngin = fopen( "tmp.png", "rb" ); im = gdImageCreateFromPng(pngin); fclose( pngin ); black = gdImageColorAllocate (im, 0, 0, 0); blue = gdImageColorAllocate (im, 0, 0, 255); red = gdImageColorAllocate (im, 224, 0, 0); // color scale ------------------------------------------------------ // 0 to 300+3*255 = 1065 for( i=0 ; i<300; i++ ) { colortempr[i] = black ; } // black to red to yellow to white // only 256 colors for image, assign 153 of them for( i=0 ; i<255; i+=5 ) { colortempr[i+300] = gdImageColorAllocate( im, i, 0, 0 ); colortempr[i+555] = gdImageColorAllocate( im, 255, i, 0 ); colortempr[i+810] = gdImageColorAllocate( im, 255, 255, i ); for( j=1 ; j<5 ; j++ ) { colortempr[i+j+300]=colortempr[i+300]; colortempr[i+j+555]=colortempr[i+555]; colortempr[i+j+810]=colortempr[i+810]; } } // right side color temperature scale for( iy=YGP_UL ; iy <= YGP_LR ; iy++ ) { itempr = (int) ( TGP_LO + ( (double) ( TGP_HI - TGP_LO ) ) * ( (double) ( iy - YGP_LR ) ) / ( (double) ( YGP_UL - YGP_LR ) ) ) ; if ( itempr < 930 ) { gdImageLine(im, XGP_LR+ 1, iy, XGP_LR+10, iy, colortempr [ itempr ] ); } } // track temperature for( ix = XGP_UL ; ix <= XGP_LR ; ix++ ) { ipoint = (int) ( ( (double) ( PPLOT - OPLOT ) ) * ( (double) ( ix - XGP_UL ) ) / ( (double) ( XGP_LR - YGP_UL ) ) ) ; // forward track // points between (-)OPLOT and PPLOT , minus the index ip = rmod( OPLOT+ipoint-index ); itempr = (int) ( rtemp[ ip ] ); if ( itempr >= MAXTEMP ) { itempr = MAXTEMP-1 ; } gdImageLine(im, ix, YGP_UL+ 7, ix, YGP_UL-10, colortempr[ itempr ] ); // reverse track // points between RL/2+PPLOT and RL/2-OPLOT, minus the index jpoint = RL/2+PPLOT-ipoint ; ip = rmod( OPLOT+jpoint-index ); itempr = (int) ( rtemp[ ip ] ); if ( itempr >= MAXTEMP ) { itempr = MAXTEMP-1 ; } gdImageLine(im, ix, YGP_UL+43, ix, YGP_UL+50, colortempr[ itempr ] ); } // turnarounds for( ipoint = PPLOT ; ipoint < ((RL/2)); ipoint++ ) { for( tsmooth=0.9 ; tsmooth > -0.01 ; tsmooth -= 0.1 ) { // smoothing tangle = PI* ( tsmooth + (double)(ipoint-PPLOT)) / ( (double)((RL/2)-PPLOT)) ; // east turnaround ip = rmod( ipoint+OPLOT-index ); x1 = XGP_LR +18*sin(tangle); y1 = YGP_UL+25 -18*cos(tangle); x2 = XGP_LR +30*sin(tangle); y2 = YGP_UL+20 -30*cos(tangle); itempr = (int) ( rtemp[ ip ] ); gdImageLine(im, x1, y1, x2, y2, colortempr[ itempr ] ); // west turnaround ip = rmod( RL/2+ipoint+OPLOT-index ); x1 = XGP_UL -18*sin(tangle); y1 = YGP_UL+25 +18*cos(tangle); x2 = XGP_UL -30*sin(tangle); y2 = YGP_UL+20 +30*cos(tangle); itempr = (int) ( rtemp[ ip ] ); gdImageLine(im, x1, y1, x2, y2, colortempr[ itempr ] ); } } // payloads on top of forward track for( i = 0 ; i < npay ; i++ ) { ix = x_to_i( xpay[i] ); payload[0].x = ix ; payload[0].y = YGP_UL - 10 ; payload[1].x = ix-15 ; payload[1].y = YGP_UL - 25 ; payload[2].x = ix+15 ; payload[2].y = YGP_UL - 25 ; gdImageFilledPolygon( im, payload, 3, blue ); } // payload heat tracks for( i = 0 ; i < PN ; i++ ) { ix = x_to_i( xht0[i] ); hx = x_to_i( xht1[i] ); if ( hx > 0 ) { if ( ix < XGP_LR ) { if ( hx > XGP_LR ) { hx = XGP_LR ; } iy = YGP_LR - (i+1)*YHSTEP ; gdImageFilledRectangle( im, ix, iy, hx, iy+YHBAND , colortempr[ cht[i] ] ); } } } // ------------------------------------------------------------------ sprintf( plotlabel, "%4d/%4d sec,%3d vehicles,%4dx animation speedup\n", ((int)time), TIME, PN, speedup ); gdImageStringFT( im, NULL, // imagespace, bounding box black , FNT, FSZ, 0.0, // color, font, fontsize, angle XGP_UL+2, YGP_UL+36, plotlabel ); // x, y, text pngout = fopen( "tmpl.png", "wb"); gdImagePngEx( im, pngout, 1 ); gdImageDestroy(im); fclose(pngout); // move temp file to directory sprintf( filename, "%sdir/b%04d.png", NAME, iplot ); rename( "tmpl.png", filename ); } // make flash video sprintf( png2swf, "png2swf -o %s.swf -r%3d %sdir/*.png", NAME, RATE, NAME ); system( png2swf ); return(0); }