00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <stdio.h>
00018 #include <stdlib.h>
00019 #include <string.h>
00020 #include <math.h>
00021 #include "defs.h"
00022 #include "athena.h"
00023 #include "prototypes.h"
00024
00025 static Real **data=NULL;
00026
00027
00028
00029
00030
00031
00032 static void compute_rgb(double data, double min, double max, int *pR, int *pG,
00033 int *pB, OutputS *p);
00034
00035
00036
00037
00038
00039 void output_ppm(MeshS *pM, OutputS *pOut)
00040 {
00041 GridS *pGrid;
00042 FILE *pfile;
00043 char *fname,*plev=NULL,*pdom=NULL;
00044 char levstr[8],domstr[8];
00045 int nl,nd,nx1,nx2,i,j;
00046 Real dmin, dmax;
00047 int red,green,blue;
00048
00049
00050 if (pOut->ndim != 2) {
00051 ath_error("[output_ppm:] Data must be 2D\n");
00052 }
00053
00054
00055
00056 for (nl=0; nl<(pM->NLevels); nl++){
00057 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00058 if (pM->Domain[nl][nd].Grid != NULL){
00059
00060
00061 if ((pOut->nlevel == -1 || pOut->nlevel == nl) &&
00062 (pOut->ndomain == -1 || pOut->ndomain == nd)){
00063 pGrid = pM->Domain[nl][nd].Grid;
00064
00065
00066
00067
00068
00069 data = OutData2(pGrid,pOut,&nx1,&nx2);
00070 if (data != NULL){
00071
00072
00073
00074
00075 if (nl>0) {
00076 plev = &levstr[0];
00077 sprintf(plev,"lev%d",nl);
00078 }
00079 if (nd>0) {
00080 pdom = &domstr[0];
00081 sprintf(pdom,"dom%d",nd);
00082 }
00083
00084 if((fname = ath_fname(plev,pM->outfilename,plev,pdom,num_digit,
00085 pOut->num,pOut->id,"ppm")) == NULL){
00086 ath_error("[output_ppm]: Error constructing filename\n");
00087 }
00088
00089
00090 if((pfile = fopen(fname,"w")) == NULL){
00091 ath_error("[output_ppm]: Unable to open ppm file %s\n",fname);
00092 }
00093
00094
00095 minmax2(data,nx2,nx1,&dmin,&dmax);
00096 pOut->gmin = MIN(dmin,pOut->gmin);
00097 pOut->gmax = MAX(dmax,pOut->gmax);
00098
00099 fprintf(pfile,"P6\n");
00100 fprintf(pfile,"# dmin= %.7e, dmax= %.7e, gmin= %.7e, gmax= %.7e\n",
00101 dmin,dmax,pOut->gmin,pOut->gmax);
00102 fprintf(pfile,"%d %d\n255\n",nx1,nx2);
00103
00104
00105 if (pOut->sdmin != 0) dmin = pOut->dmin;
00106 if (pOut->sdmax != 0) dmax = pOut->dmax;
00107
00108 for (j=nx2-1; j>=0; j--) {
00109 for (i=0; i<nx1; i++) {
00110 compute_rgb(data[j][i],dmin,dmax,&red,&green,&blue,pOut);
00111 fputc(red,pfile);
00112 fputc(green,pfile);
00113 fputc(blue,pfile);
00114 }
00115 }
00116
00117
00118 fclose(pfile);
00119 free_2d_array(data);
00120 free(fname);
00121 data = NULL;
00122 }
00123 }}
00124 }
00125 }
00126 }
00127
00128
00129
00130
00131
00132
00133
00134 static void compute_rgb(double data, double min, double max,
00135 int *R, int *G, int *B, OutputS *pOut)
00136 {
00137 int i;
00138 float x, *rgb = pOut->rgb, *der = pOut->der;
00139
00140 if (rgb == 0) {
00141 *R = *G = *B = (data > max ? 255 : 0);
00142 return;
00143 }
00144
00145 if (min==max) {
00146 *R = *G = *B = (data > max ? 255 : 0);
00147 return;
00148 }
00149 #if 1
00150 x = (data-min)*255.0/(max-min);
00151 if (x<=0.0 || x>=255.0) {
00152 i= (x <= 0.0 ? 0 : 255);
00153 *R = (int) (rgb[i*3+0] * 255.0);
00154 *G = (int) (rgb[i*3+1] * 255.0);
00155 *B = (int) (rgb[i*3+2] * 255.0);
00156 return;
00157 }
00158 i = (int) x;
00159 *R = (int) ((rgb[i*3+0] + (x-i)*der[i*3+0])*255.0);
00160 *G = (int) ((rgb[i*3+1] + (x-i)*der[i*3+1])*255.0);
00161 *B = (int) ((rgb[i*3+2] + (x-i)*der[i*3+2])*255.0);
00162 #else
00163 i = (int) ((data-min)*255.0/(max-min));
00164 if (i<0) i=0;
00165 if (i>255) i=255;
00166 *R = (int) (rgb[i*3+0] * 255.0);
00167 *G = (int) (rgb[i*3+1] * 255.0);
00168 *B = (int) (rgb[i*3+2] * 255.0);
00169 #endif
00170 }