diff --git a/HEN_HOUSE/estar/estarMainCalc.cpp b/HEN_HOUSE/estar/estarMainCalc.cpp
index 58cba5e4d..d99755676 100644
--- a/HEN_HOUSE/estar/estarMainCalc.cpp
+++ b/HEN_HOUSE/estar/estarMainCalc.cpp
@@ -13,13 +13,83 @@
 #include "modules/solverHelpers/cutoff.cpp"
 using namespace std;
 
+string getFileNameWithoutExtension(const string &s) {
+    char sep = '/';
+#ifdef _WIN32
+    sep = '\\';
+#endif
+
+    size_t i = s.rfind(sep, s.length());
+    if (i != string::npos) {
+        string filename = s.substr(i+1, s.length() - i);
+        size_t lastindex = filename.find_last_of(".");
+        string rawname = filename.substr(0, lastindex);
+        return (rawname);
+    }
+    return ("");
+}
+
+// Output a density correction file
+void outputDensityFile(float mediaDensity, double *densityCorr, double *enGrid, float *meanIval, formula_calc fc, string outputFilename) {
+    // Remove spaces from the end because it was a fixed length fortran string
+    outputFilename = outputFilename.erase(outputFilename.find_last_not_of(" \n\r\t")+1);
+    if (outputFilename.empty()) {
+        return;
+    }
+
+    cout << "Writing density correction file: '" << outputFilename.c_str() << "'" << endl;
+
+    std::ofstream f(outputFilename.c_str());
+
+    f << getFileNameWithoutExtension(outputFilename).c_str() << endl;
+    f << setprecision(8);
+    f << "113 " << *meanIval << " " << mediaDensity << " " << fc.mmax << endl;
+
+    // Output the atomic numbers and mass fractions
+    unsigned int k = 0;
+    while (k < fc.mmax) {
+        f << fc.jz[k] << " " << fc.wt[k];
+
+        // Put a new line every 6 elements, or after the last element
+        if ((k+1) % 6 == 0 || k+1 == fc.mmax) {
+            f << endl;
+        }
+        else {
+            f << " ";
+        }
+
+        k += 1;
+    }
+
+    // Set formatting for the density data
+    f << scientific << showpoint;
+
+    // Output the energy grid and density effect corrections
+    for (int i = 0; i < 113; i++) {
+        f << setprecision(2);
+        f << enGrid[i];
+        f << setprecision(3);
+        f << "," << densityCorr[i];
+
+        // Put a new line every 4 values, or after the last value
+        if ((i+1) % 4 == 0 || i+1 == 113) {
+            f << endl;
+        }
+        else {
+            f << "  ";
+        }
+    }
+
+    f.close();
+}
+
 /*
     This function computes the density correction factors by using the processed input arrays and variables
     from estarCalc.cpp. The density correction factors are stored in the variable densityCorr.
 */
 
 int estarCalculation(int isCompound, int NEP, float mediaDensity, string *elementArray, double *massFraction,
-                     float *numOfAtoms, double *densityCorr, double *enGrid, float *meanIval, float *ipotval, int mediaNum) {
+                     float *numOfAtoms, double *densityCorr, double *enGrid, float *meanIval, float *ipotval, int mediaNum, string outputFilename) {
     //------------------------------------------------//
     int knmat;
 
@@ -376,6 +446,10 @@ int estarCalculation(int isCompound, int NEP, float mediaDensity, string *elemen
     cout << "\n";
     cout << "Density correction factors calculated by ESTAR for medium " << mediaNum << ".\n";
     cout << "-------------------------\n";
+
+    // output a density correction file
+    outputDensityFile(mediaDensity, densityCorr, enGrid, meanIval, fc, outputFilename);
+
     return 0;
 
 };
diff --git a/HEN_HOUSE/estar/estarPreProcess.cpp b/HEN_HOUSE/estar/estarPreProcess.cpp
index ba6f19b39..b3f360407 100644
--- a/HEN_HOUSE/estar/estarPreProcess.cpp
+++ b/HEN_HOUSE/estar/estarPreProcess.cpp
@@ -123,7 +123,9 @@ extern "C" int estar_(char *formulaStr,
                       int *ISCOMP,
                       float *meanIval,
                       float *ipotval,
-                      int *mediaID) {
+                      int *mediaID,
+                      char *outputFilename
+                     ) {
     string mainFormula;
     string mainFormula_temp_1;
     string mainFormula_temp_2;
@@ -148,7 +150,7 @@ extern "C" int estar_(char *formulaStr,
     double mediumDensity = *mediaDensity;
     // call the main estar calculation files.
     estarCalculation(isCompInt, nepInt, mediumDensity, estarFormulaArrayInput, estarWeightArrayInput,
-                     numOfAtoms, densityCorr, enGrid, meanIval, ipotval, mediaNum);
+                     numOfAtoms, densityCorr, enGrid, meanIval, ipotval, mediaNum, string(outputFilename));
 
     return 0;
 }
diff --git a/HEN_HOUSE/src/get_media_inputs.mortran b/HEN_HOUSE/src/get_media_inputs.mortran
index 2fbcefa2a..5c92f32be 100644
--- a/HEN_HOUSE/src/get_media_inputs.mortran
+++ b/HEN_HOUSE/src/get_media_inputs.mortran
@@ -668,13 +668,13 @@ COMIN/MEDINP,ELECIN,THRESH,ELEMTB,USEFUL/;
 $INTEGER  ival,ival_media,ival_medfile,i,j,k,ival_ae,ival_ue,ival_ap,ival_up,
           ival_rho,ival_elements,ival_rhoz,ival_iunrst,ival_iaprim,ival_gasp,
           ival_pz,ival_sterncid,ival_ipot,ival_localdensity,
-          ival_densityfile,medfile_error,ival_outfile,
+          ival_densityfile,medfile_error,ival_outfile,ival_densityoutputfile,
           egs_open_file,lnblnk1,i_medfile,egs_get_unit,i_mederr,mindex,eindex,
           i_density,i01,length,i_outfile;
 $REAL ecut_min, pcut_min;
 
 $LOGICAL  medfile_specified,densityfile_specified,elements_specified,
-          outfile_specified($MXMED);
+          outfile_specified($MXMED), densityoutputfile_specified;
 $LOGICAL  iunrst_specified,stern_specified,iaprim_specified,
           gasp_specified,rho_specified,start_delim_found,end_delim_found,
           spec_by_pz,spec_by_rhoz,df_if_elem_mismatch($MXMED),
@@ -701,7 +701,7 @@ DATA ETAB/1.,1.25,1.5,1.75,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,7.,8.,9./;
 
 character*24 medium_name,med_tmp,sterncid_tmp;
 character*256 density_file,material_file,tmp_string,
-              spoutput_file($MXMED);
+              spoutput_file($MXMED), density_outputfile;
 character*80 text_string, text_save, title;
 character*80 delim_start,delim_end;
 character*1 blank;
@@ -1079,6 +1079,19 @@ DO i=1,NMED[
       localdensity_tmp=rho_tmp;
    ]
 
+   ival=ival+1;
+   ival_densityoutputfile=ival;
+   values_sought(ival)='output density file';
+   type(ival)          = 2;
+   nvalue(ival)=1;
+   nmin=ival_densityoutputfile;nmax=ival_densityoutputfile;
+   CALL GET_INPUT;
+   IF(error_flags(ival_densityoutputfile)=0)[
+       density_outputfile=char_value(ival_densityoutputfile,1);
+       densityoutputfile_specified=.true.;
+       $WRITE_MEDERR(' Density correction file will be output.');
+   ]
+
    "now see what else is specified in the .egsinp file"
    ival=ival+1;
    ival_ipot=ival;
@@ -1662,6 +1675,12 @@ DO i=1,NMED[
        inpgasp(i)=gasp_tmp;
        inpdensity_file(i)=density_file;
 
+       "Only output the density correction file for the first job if parallel"
+       IF(densityoutputfile_specified &
+            (n_parallel=0 | i_parallel=first_parallel)) [
+         outdensity_file(i)=density_outputfile;
+       ]
+
     ]
     ELSE[
       $WRITE_MEDERR(' Error: Medium ',medium_name,' not correctly defined.');
@@ -1679,6 +1698,7 @@ IF ( ounit <= 0 ) return;
 
 IF(is_pegsless)[
 
+
 "show common data"
 
 write(ounit,*);
@@ -1703,13 +1723,13 @@ DO i=1,nmed[
   write(ounit,'(a,24a1)')'   Medium: ',(media(j,i),j=1,24);
   write(ounit,'(a,24a1)')' Sterncid: ',(inpstrn(j,i),j=1,24);
   write(ounit,'(a,1p,e14.5,a)')' rho (bulk density): ',rho(i),' g/cm^3';
-  IF(~densityfile_specified)[
+  IF(epstfl(i)~=1)[
  write(ounit,'(a,1p,e14.5,a)')'      local density: ',localDensity(i),' g/cm^3';
   ] ELSE [
     write(ounit,'(a)')'      local density: using .density file';
   ]
   IF(~(ipot(i)=-1))[
-    write(ounit,'(a,1p,e14.5,a)')'     I-value: ',ipot(i),' MeV';
+    write(ounit,'(a,1p,e14.5,a)')' I-value: ',ipot(i),' MeV';
   ]
   write(ounit,'(a,24a4)')' Elements: ',(inpasym(i,j),j=1,nne(i));
   write(ounit,'(a,1p,12e14.5)')'    rhoz: ',(rhoz(i,j),j=1,nne(i));
@@ -1717,6 +1737,10 @@ DO i=1,nmed[
   write(ounit,'(a,i5)')'  iunrst: ',iunrst(i);
   write(ounit,'(a,i5)')'  iaprim: ',iaprim(i);
   write(ounit,'(a,1p,e14.5,a)')'    gasp: ',inpgasp(i),' atm.';
+  IF(epstfl(i)~=1)[
+    write(ounit,*)' output density correction file: ',
+        $cstring(outdensity_file(i));
+  ]
   IF(epstfl(i)=1)[
     write(ounit,*)' density correction file: ',
         $cstring(inpdensity_file(i));
diff --git a/HEN_HOUSE/src/pegs4_macros.mortran b/HEN_HOUSE/src/pegs4_macros.mortran
index abe99ed2a..351a35c30 100644
--- a/HEN_HOUSE/src/pegs4_macros.mortran
+++ b/HEN_HOUSE/src/pegs4_macros.mortran
@@ -168,7 +168,8 @@ REPLACE {;COMIN/ELEMTB/;} WITH
 "1 means the medium is a compound."
 REPLACE {;COMIN/MIXDAT/;} WITH
    {;COMMON/MIXDAT/NEP,LMED,PZP($MXEL),ZELEMP($MXEL),WAP($MXEL),RHOZP($MXEL),
-     GASPP,EZ,TPZ,IDSTRN(24),ELEM_NAMES($MXEL),ISCOMP,IPOTVAL,LDENSITY,MEDID;
+     GASPP,EZ,TPZ,IDSTRN(24),ELEM_NAMES($MXEL),ISCOMP,IPOTVAL,LDENSITY,MEDID,
+     OUTDENSITY;
      $INTEGER NEP,LMED;
      $REAL4 PZP,ZELEMP,WAP,RHOZP,GASPP,EZ,TPZ;
      $TYPE IDSTRN;
@@ -176,6 +177,7 @@ REPLACE {;COMIN/MIXDAT/;} WITH
      $INTEGER ISCOMP;
      $REAL4 IPOTVAL, LDENSITY;
      $INTEGER MEDID;
+     character*256 OUTDENSITY;
    }
 
 REPLACE {;COMIN/DBRPR/;} WITH
@@ -192,8 +194,9 @@ REPLACE {;COMIN/MEDINP/;} WITH
    {;COMMON/MEDINP/inpdensity_file($MXMED),inpasym($MXMED,$MXEL),
      inpstrn(24,$MXMED),pz4($MXMED,$MXEL),
      rhoz4($MXMED,$MXEL),wa4($MXMED,$MXEL),inpgasp($MXMED),
-     iscompound($MXMED),ipot($MXMED),localDensity($MXMED);
-     character*256 inpdensity_file;
+     iscompound($MXMED),ipot($MXMED),localDensity($MXMED),
+     outdensity_file($MXMED);
+     character*256 inpdensity_file,outdensity_file;
      $TYPE inpasym,inpstrn;
      $REAL4 pz4,rhoz4,wa4,inpgasp;
      $INTEGER iscompound;
@@ -235,6 +238,7 @@ DO IM=1,NMED[
   ISCOMP = iscompound(IM);
   IPOTVAL = ipot(IM);
   LDENSITY = localDensity(IM);
+  OUTDENSITY = outdensity_file(IM);
   IUNRSTP=IUNRST(IM);
   IAPRIMP=IAPRIM(IM);
   EPSTFLP=EPSTFL(IM);
diff --git a/HEN_HOUSE/src/pegs4_routines.mortran b/HEN_HOUSE/src/pegs4_routines.mortran
index 02a890a7c..21f0bb0e8 100644
--- a/HEN_HOUSE/src/pegs4_routines.mortran
+++ b/HEN_HOUSE/src/pegs4_routines.mortran
@@ -1705,7 +1705,7 @@ IPOTVAL = IPOTVAL * 1000000;
 "Here we call estar"
 "The parameters are described in estarCalc.cpp"
 CALL estar(ELEM_NAMES,RHOZP,PZP,LDENSITY,densityCorr,
-enGrid,NEP,ISCOMP,IEV,IPOTVAL,MEDID);
+enGrid,NEP,ISCOMP,IEV,IPOTVAL,MEDID,OUTDENSITY);
 
 "After the call to ESTAR, the remainder"
 "of the calculation is similiar to the EPSTFL=1"