Skip to content

Commit

Permalink
Add option to output estar density file
Browse files Browse the repository at this point in the history
  • Loading branch information
@rtownson
rtownson committed Jul 25, 2022
1 parent 953630e commit 914a03a
Show file tree
Hide file tree
Showing 5 changed files with 116 additions and 12 deletions.
76 changes: 75 additions & 1 deletion HEN_HOUSE/estar/estarMainCalc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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;

};
6 changes: 4 additions & 2 deletions HEN_HOUSE/estar/estarPreProcess.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
}
Expand Down
34 changes: 29 additions & 5 deletions HEN_HOUSE/src/get_media_inputs.mortran
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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.');
Expand All @@ -1679,6 +1698,7 @@ IF ( ounit <= 0 ) return;

IF(is_pegsless)[


"show common data"

write(ounit,*);
Expand All @@ -1703,20 +1723,24 @@ 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));
write(ounit,'(a,1p,12e14.5)')' pz: ',(pz(i,j),j=1,nne(i));
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));
Expand Down
10 changes: 7 additions & 3 deletions HEN_HOUSE/src/pegs4_macros.mortran
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,14 +168,16 @@ 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;
CHARACTER*2 ELEM_NAMES;
$INTEGER ISCOMP;
$REAL4 IPOTVAL, LDENSITY;
$INTEGER MEDID;
character*256 OUTDENSITY;
}

REPLACE {;COMIN/DBRPR/;} WITH
Expand All @@ -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;
Expand Down Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion HEN_HOUSE/src/pegs4_routines.mortran
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down

0 comments on commit 914a03a

Please sign in to comment.