diff --git a/residue.cpp b/residue.cpp
index 35833fc..de48fde 100644
--- a/residue.cpp
+++ b/residue.cpp
@@ -1054,13 +1054,15 @@ void printSpaceTimeHeader(UserData data)
 		fprintf((data->output), "%15s\t",data->gas->speciesName(k-1).c_str());
 	}
 	fprintf((data->output), "%15s\t","Pressure(Pa)");
-	fprintf((data->output), "%15s\n","HRR(J)");
+	fprintf((data->output), "%15s\t","HRR(J)");
+	fprintf((data->output), "%15s\n","rho(kg/m3)");
 
 }
 
 void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
 {
 	double *ydata,*psidata;
+	double rho;
   	ydata    = N_VGetArrayPointer_OpenMP(*y);
 
 	if(data->adaptiveGrid){
@@ -1070,6 +1072,8 @@ void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
 	}
 
 	for (size_t i = 0; i < data->npts; i++) {
+		setGas(data,ydata,i+1,GAS) ; 
+		rho = data->gas->density() ;
 		fprintf(output, "%15.6e\t%15.6e\t",psi(i+1),t);
 		if(i==0){
 			fprintf(output, "%15.6e\t",psi(2)-psi(1));
@@ -1081,6 +1085,7 @@ void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
 			fprintf(output, "%15.9e\t",ydata[j+i*data->nvar]);
 		}
 		fprintf(output,"%15.6e\t",data->HRRdata[i]);
+		fprintf(output,"%15.9e\t",rho);
 		fprintf(output, "\n");
 	}
 	fprintf(output, "\n\n");