import rock.*;

import java.awt.*;
import java.io.*;

import javax.swing.*;

public class MeltsOutput {
  private rock.dRock rockRef;
  
  public PhaseHandler(rock.dRock rockRef) {
    this.rockRef = rockRef;
  }

// output -> melts.out
// liquidFile -> melts-liquid.tbl

    fprintf(tableLiq, "Index,T (C),P (kbars),log(10) f O2");
    fprintf(tableLiq, ",liq mass (gm),liq rho (gm/cc)");
    for (i=0; i<nc; i++) fprintf(tableLiq, ",wt%% %s", bulkSystem[i].label);
    fprintf(tableLiq, ",liq G (kJ),liq H (kJ),liq S (J/K),liq V (cc),liq Cp (J/K)");
    for (i=0; i<nlc; i++) fprintf(tableLiq, ",activity %s", liquid[i].label);
    fprintf(tableLiq, ",liq vis (log 10 poise),sol mass (gm),sol rho (gm/cc)");
    fprintf(tableLiq, ",sol G (kJ),sol H (kJ),sol S (J/K),sol V (cc),sol Cp (J/K)");
    fprintf(tableLiq, "\n");

  fprintf(output, "\n**********----------**********\nTitle: %s\n\n", "Dummy Title");
  fprintf(output, "T = %.2f (C)  P = %.3f (kbars)  log(10) f O2 = %.2f  ",
    silminState->T-273.15, silminState->P/1000.0, silminState->fo2);
  fprintf(output, "delta HM = %.2f  NNO = %.2f  QFM = %.2f  COH = %.2f  IW = %.2f\n\n",
    silminState->fo2 - getlog10fo2(silminState->T, silminState->P, FO2_HM), 
    silminState->fo2 - getlog10fo2(silminState->T, silminState->P, FO2_NNO), 
    silminState->fo2 - getlog10fo2(silminState->T, silminState->P, FO2_QFM),
    silminState->fo2 - getlog10fo2(silminState->T, silminState->P, FO2_COH), 
    silminState->fo2 - getlog10fo2(silminState->T, silminState->P, FO2_IW));
  
  fprintf(output, "Constraint Flags: ");
  
  if      (silminState->fo2Path == FO2_HM)     fprintf(output, "fO2 path = HM  ");
  else if (silminState->fo2Path == FO2_QFM)    fprintf(output, "fO2 path = QFM  ");
  else if (silminState->fo2Path == FO2_COH)    fprintf(output, "fO2 path = C-COH  ");
  else if (silminState->fo2Path == FO2_NNO)    fprintf(output, "fO2 path = NNO  ");
  else if (silminState->fo2Path == FO2_IW)     fprintf(output, "fO2 path = IW  ");
  else if (silminState->fo2Path == FO2_QFM_P3) fprintf(output, "fO2 path = QFM+3  ");
  else if (silminState->fo2Path == FO2_QFM_P2) fprintf(output, "fO2 path = QFM+2  ");
  else if (silminState->fo2Path == FO2_QFM_P1) fprintf(output, "fO2 path = QFM+1  ");
  else if (silminState->fo2Path == FO2_QFM_M1) fprintf(output, "fO2 path = QFM-1  ");
  else if (silminState->fo2Path == FO2_QFM_M2) fprintf(output, "fO2 path = QFM-2  ");
  else if (silminState->fo2Path == FO2_QFM_M3) fprintf(output, "fO2 path = QFM-3  ");
  else if (silminState->fo2Path == FO2_QFM_M4) fprintf(output, "fO2 path = QFM-4  ");
  else if (silminState->fo2Path == FO2_QFM_M5) fprintf(output, "fO2 path = QFM-5  ");
  else if (silminState->fo2Path == FO2_QFM_M6) fprintf(output, "fO2 path = QFM-6  ");
  else if (silminState->fo2Path == FO2_QFM_M7) fprintf(output, "fO2 path = QFM-7  ");
  else if (silminState->fo2Path == FO2_QFM_M8) fprintf(output, "fO2 path = QFM-8  ");
  else if (silminState->fo2Path == FO2_QFM_M9) fprintf(output, "fO2 path = QFM-9  ");
  
  if (silminState->fractionateSol) fprintf(output, "Fractionate Solids  ");
  if (silminState->fractionateLiq) fprintf(output, "Fractionate Liquids  ");
  if (silminState->multipleLiqs)   fprintf(output, "Multiple Liquids  ");
  if (silminState->isenthalpic)    fprintf(output, "Isenthalpic Path  ");
  if (silminState->isentropic)     fprintf(output, "Isentropic Path  ");
  if (silminState->isochoric)      fprintf(output, "Isochoric Path  ");
  if (silminState->assimilate)     fprintf(output, "Assimilation  ");

  fprintf(output, "\n\n");

  if (hasLiquid) {
    gLiq = 0.0; hLiq = 0.0; sLiq = 0.0; vLiq = 0.0; cpLiq = 0.0; mLiq = 0.0;
    for (nl=0; nl<silminState->nLiquidCoexist; nl++) {
      conLiq(SECOND, THIRD, silminState->T, silminState->P, NULL, silminState->liquidComp[nl], r, NULL, NULL, NULL, NULL);

      gmixLiq (FIRST, silminState->T, silminState->P, r, &gibbsEnergy,  NULL, NULL);
      hmixLiq (FIRST, silminState->T, silminState->P, r, &enthalpy,     NULL);
      smixLiq (FIRST, silminState->T, silminState->P, r, &entropy,      NULL, NULL, NULL);
      vmixLiq (FIRST, silminState->T, silminState->P, r, &volume,       NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
      cpmixLiq(FIRST, silminState->T, silminState->P, r, &heatCapacity, NULL, NULL);
      visLiq  (FIRST, silminState->T, silminState->P, r, &viscosity);

      for (i=0, moles=0.0; i<nlc; i++) moles +=  (silminState->liquidComp)[nl][i];
      gibbsEnergy *= moles; enthalpy	 *= moles; entropy *= moles;
      volume	  *= moles; heatCapacity *= moles;

      for (i=0; i<nlc; i++) {
    	gibbsEnergy  += (silminState->liquidComp)[nl][i]*(liquid[i].cur).g;
    	enthalpy     += (silminState->liquidComp)[nl][i]*(liquid[i].cur).h;
    	entropy      += (silminState->liquidComp)[nl][i]*(liquid[i].cur).s;
    	volume       += (silminState->liquidComp)[nl][i]*(liquid[i].cur).v;
    	heatCapacity += (silminState->liquidComp)[nl][i]*(liquid[i].cur).cp;
      }

      for (i=0, oxSum=0.0; i<nc; i++) {
    	for (j=0, oxVal[i]=0.0; j<nlc; j++) oxVal[i] += (liquid[j].liqToOx)[i]*(silminState->liquidComp)[nl][j];
    	oxVal[i] *= bulkSystem[i].mw;
    	oxSum	 += oxVal[i];
      }
      if (oxSum != 0.0) for (i=0; i<nc; i++) oxVal[i] /= oxSum;

      gLiq += gibbsEnergy; hLiq += enthalpy; sLiq += entropy; vLiq += volume; cpLiq += heatCapacity; mLiq += oxSum;

      fprintf(output, "%-15.15s  mass = %.2f (gm)  density = %.2f (gm/cc)", "Liquid", oxSum, (volume == 0.0) ? 0.0 : oxSum/(10.0*volume));
      fprintf(output, "  viscosity = %.2f (log 10 poise)", viscosity);
      fprintf(output, "     (analysis in wt %%)\n");

      fprintf(output, " 		");
      fprintf(output, "G = %.2f (J)  ",    gibbsEnergy);
      fprintf(output, "H = %.2f (J)  ",    enthalpy);
      fprintf(output, "S = %.2f (J/K)  ",  entropy);
      fprintf(output, "V = %.2f (cc)  ",   volume*10.0);
      fprintf(output, "Cp = %.2f (J/K)  ", heatCapacity);
      fprintf(output, "\n");
 
      fprintf(output, "     ");
      for (i=0; i<nc; i++) fprintf(output, "%7.7s", bulkSystem[i].label);
      fprintf(output, "\n");
      fprintf(output, "     ");
      for (i=0; i<nc; i++) fprintf(output, "  %5.2f", oxVal[i]*100.0);
      fprintf(output, "\n");

      fprintf(tableLiq, "%d,%.2f,%.3f,%.3f", rowIndex, silminState->T-273.15, silminState->P/1000.0, silminState->fo2);
      fprintf(tableLiq, ",%.4f,%.4f", oxSum, (volume == 0.0) ? 0.0 : oxSum/(10.0*volume));
      for (i=0; i<nc; i++) fprintf(tableLiq, ",%.4f", oxVal[i]*100.0);
      fprintf(tableLiq, ",%.3f,%.3f,%.3f,%.3f,%.3f", gibbsEnergy/1000.0, enthalpy/1000.0, entropy, volume*10.0, heatCapacity);
      actLiq(FIRST, silminState->T, silminState->P, r, m, NULL, NULL, NULL);
      for (i=0; i<nlc; i++) fprintf(tableLiq, ",%.7e", m[i]);
      fprintf(tableLiq, ",%.4f", viscosity);
  
    }
  } else {
    gLiq = 0.0; hLiq = 0.0; sLiq = 0.0; vLiq = 0.0; cpLiq = 0.0; mLiq = 0.0;
  }

  for (j=0, totalMass=0.0, totalGibbsEnergy=0.0, totalEnthalpy=0.0,
    totalEntropy=0.0, totalVolume=0.0, totalHeatCapacity=0.0; j<npc; j++) 
  for (ns=0; ns<(silminState->nSolidCoexist)[j]; ns++) {
    if (solids[j].na == 1) {
      mass               = (silminState->solidComp)[j][ns]*solids[j].mw;
      gibbsEnergy        = (silminState->solidComp)[j][ns]*(solids[j].cur).g;
      enthalpy           = (silminState->solidComp)[j][ns]*(solids[j].cur).h;
      entropy            = (silminState->solidComp)[j][ns]*(solids[j].cur).s;
      volume             = (silminState->solidComp)[j][ns]*(solids[j].cur).v;
      heatCapacity       = (silminState->solidComp)[j][ns]*(solids[j].cur).cp;
      totalMass         += (silminState->solidComp)[j][ns]*solids[j].mw;
      totalGibbsEnergy  += (silminState->solidComp)[j][ns]*(solids[j].cur).g;
      totalEnthalpy     += (silminState->solidComp)[j][ns]*(solids[j].cur).h;
      totalEntropy      += (silminState->solidComp)[j][ns]*(solids[j].cur).s;
      totalVolume       += (silminState->solidComp)[j][ns]*(solids[j].cur).v;
      totalHeatCapacity += (silminState->solidComp)[j][ns]*(solids[j].cur).cp;
      fprintf(output, "\n%-15.15s  mass = %.2f (gm)  density = %.2f (gm/cc)", solids[j].label, mass, (volume == 0.0) ? 0.0 : mass/(10.0*volume));
      fprintf(output, "\n");
      fprintf(output, "                 %s\n", solids[j].formula); 
      fprintf(output, "                 "); 
      fprintf(output, "G = %.2f (J)  ", gibbsEnergy);
      fprintf(output, "H = %.2f (J)  ", enthalpy);
      fprintf(output, "S = %.2f (J/K)  ", entropy);
      fprintf(output, "V = %.2f (cc)  ", volume*10.0);
      fprintf(output, "Cp = %.2f (J/K)  ", heatCapacity);
      fprintf(output, "\n");
      
      // make a phase.tbl file -> tableSol
      // top labels
        fprintf(tableSol[j], "Index,T (C),P (kbars),log(10) f O2");
        fprintf(tableSol[j], ",mass (gm),rho (gm/cc)");
        for (i=0; i<nc; i++) fprintf(tableSol[j], ",wt%% %s", bulkSystem[i].label);
        fprintf(tableSol[j], ",G (kJ),H (kJ),S (J/K),V (cc),Cp (J/K)");
        fprintf(tableSol[j], "\n");

      fprintf(tableSol[j], "%d,%.2f,%.3f,%.3f", rowIndex, silminState->T-273.15, silminState->P/1000.0, silminState->fo2);
      fprintf(tableSol[j], ",%.4f,%.4f", mass, (volume == 0.0) ? 0.0 : mass/(10.0*volume));
      
      for (i=0, oxSum=0.0; i<nc; i++) {
        oxVal[i]  = (solids[j].solToOx)[i]*bulkSystem[i].mw;
    	oxSum	 += oxVal[i];
      }
      if (oxSum != 0.0) for (i=0; i<nc; i++) oxVal[i] /= oxSum;
      for (i=0; i<nc; i++) fprintf(tableSol[j], ",%.4f", oxVal[i]*100.0);

      fprintf(tableSol[j], ",%.3f,%.3f,%.3f,%.3f,%.3f", gibbsEnergy/1000.0, enthalpy/1000.0, entropy, volume*10.0, heatCapacity);
      fprintf(tableSol[j], "\n");

    } else {
      char *formula;
      for (i=0, mass=0.0; i<solids[j].na; i++) {
        m[i] = (silminState->solidComp)[j+1+i][ns];
        mass += (silminState->solidComp)[j+1+i][ns]*solids[j+1+i].mw;
      }
      (*solids[j].convert)(SECOND, THIRD, silminState->T, silminState->P, NULL, m, r, NULL, NULL, NULL, NULL, NULL);
      (*solids[j].display)(FIRST, silminState->T, silminState->P, r, &formula);
      (*solids[j].gmix) (FIRST, silminState->T, silminState->P, r, &gibbsEnergy,  NULL, NULL, NULL);
      (*solids[j].hmix) (FIRST, silminState->T, silminState->P, r, &enthalpy);
      (*solids[j].smix) (FIRST, silminState->T, silminState->P, r, &entropy,      NULL, NULL);
      (*solids[j].vmix) (FIRST, silminState->T, silminState->P, r, &volume,       NULL, NULL, NULL, NULL, NULL, NULL,  NULL, NULL, NULL);
      (*solids[j].cpmix)(FIRST, silminState->T, silminState->P, r, &heatCapacity, NULL, NULL);
      gibbsEnergy  *= (silminState->solidComp)[j][ns]; 
      enthalpy     *= (silminState->solidComp)[j][ns]; 
      entropy      *= (silminState->solidComp)[j][ns];
      volume       *= (silminState->solidComp)[j][ns];
      heatCapacity *= (silminState->solidComp)[j][ns];
      for (i=0; i<solids[j].na; i++) {
        gibbsEnergy  += m[i]*(solids[j+1+i].cur).g;
        enthalpy     += m[i]*(solids[j+1+i].cur).h;
        entropy      += m[i]*(solids[j+1+i].cur).s;
        volume       += m[i]*(solids[j+1+i].cur).v;
        heatCapacity += m[i]*(solids[j+1+i].cur).cp;
      }
      totalMass         += mass;
      totalGibbsEnergy  += gibbsEnergy;
      totalEnthalpy     += enthalpy;
      totalEntropy      += entropy;
      totalVolume       += volume;
      totalHeatCapacity += heatCapacity;
      fprintf(output, "\n%-15.15s  mass = %.2f (gm)  density = %.2f (gm/cc)", solids[j].label, mass, (volume == 0.0) ? 0.0 : mass/(10.0*volume));
      fprintf(output, "     (analysis in mole %%)\n");
      fprintf(output, "                 %s\n", formula); 
      fprintf(output, "                 "); 
      fprintf(output, "G = %.2f (J)  ", gibbsEnergy);
      fprintf(output, "H = %.2f (J)  ", enthalpy);
      fprintf(output, "S = %.2f (J/K)  ", entropy);
      fprintf(output, "V = %.2f (cc)  ", volume*10.0);
      fprintf(output, "Cp = %.2f (J/K)  ", heatCapacity);
      fprintf(output, "\n");
      fprintf(output, "     "); 
      for (i=0; i<solids[j].na; i++) fprintf(output, " %13.13s", solids[j+1+i].label);
      fprintf(output, "\n");
      fprintf(output, "     "); 
      for (i=0; i<solids[j].na; i++) fprintf(output, " %13.2f", 100.0*m[i]/(silminState->solidComp)[j][ns]);
      fprintf(output, "\n");


      // make a phase.tbl file -> tableSol
      // top labels
        fprintf(tableSol[j], "Index,T (C),P (kbars),log(10) f O2");
        fprintf(tableSol[j], ",mass (gm),rho (gm/cc)");
        for (i=0; i<nc; i++) fprintf(tableSol[j], ",wt%% %s", bulkSystem[i].label);
        fprintf(tableSol[j], ",G (kJ),H (kJ),S (J/K),V (cc),Cp (J/K)");
        for (i=0; i<solids[j].na; i++) fprintf(tableSol[j], ",%13.13s", solids[j+1+i].label);
        fprintf(tableSol[j], "\n");

      fprintf(tableSol[j], "%d,%.2f,%.3f,%.3f", rowIndex, silminState->T-273.15, silminState->P/1000.0, silminState->fo2);
      fprintf(tableSol[j], ",%.4f,%.4f", mass, (volume == 0.0) ? 0.0 : mass/(10.0*volume));

      for (i=0, oxSum=0.0; i<nc; i++) {
        int k;
	for (k=0, oxVal[i]=0.0; k<solids[j].na; k++) oxVal[i] += (solids[j+1+k].solToOx)[i]*m[k]*bulkSystem[i].mw;
    	oxSum += oxVal[i];
      }
      if (oxSum != 0.0) for (i=0; i<nc; i++) oxVal[i] /= oxSum;
      for (i=0; i<nc; i++) fprintf(tableSol[j], ",%.4f", oxVal[i]*100.0);

      fprintf(tableSol[j], ",%.3f,%.3f,%.3f,%.3f,%.3f", gibbsEnergy/1000.0, enthalpy/1000.0, entropy, volume*10.0, heatCapacity);
      for (i=0; i<solids[j].na; i++) fprintf(tableSol[j], ",%.6f", m[i]/(silminState->solidComp)[j][ns]);
      fprintf(tableSol[j], "\n");
    }
  }

  fprintf(output, "\n%-15.15s  mass = %.2f (gm)  density = %.2f (gm/cc)\n", "Total solids", totalMass, (totalVolume == 0.0) ? 0.0 : totalMass/(10.0*totalVolume));
  fprintf(output, "                 "); 
  fprintf(output, "G = %.2f (J)  ", totalGibbsEnergy);
  fprintf(output, "H = %.2f (J)  ", totalEnthalpy);
  fprintf(output, "S = %.2f (J/K)  ", totalEntropy);
  fprintf(output, "V = %.2f (cc)  ", totalVolume*10.0);
  fprintf(output, "Cp = %.2f (J/K)  ", totalHeatCapacity);
  fprintf(output, "\n");

  if (hasLiquid) {
    fprintf(tableLiq, ",%.4f,%.4f", totalMass, (totalVolume == 0.0) ? 0.0 : totalMass/(10.0*totalVolume));
    fprintf(tableLiq, ",%.3f,%.3f,%.3f,%.3f,%.3f", totalGibbsEnergy/1000.0, totalEnthalpy/1000.0, totalEntropy, totalVolume*10.0, totalHeatCapacity);
    fprintf(tableLiq, "\n");
  }

  if (silminState->isenthalpic && (silminState->refEnthalpy == 0.0)) silminState->refEnthalpy = hLiq+totalEnthalpy;
  if (silminState->isentropic  && (silminState->refEntropy  == 0.0)) silminState->refEntropy  = sLiq+totalEntropy;
  if (silminState->isochoric   && (silminState->refVolume   == 0.0)) silminState->refVolume   = vLiq+totalVolume;

  if (silminState->fractionateSol || silminState->fractionateLiq) 
    fprintf(output, "\nSummary of all fractionated phases: (total mass = %.2f grams)\n", silminState->fracMass);

  if (silminState->fractionateSol) {
    for (j=0; j<npc; j++) for (ns=0; ns<(silminState->nFracCoexist)[j]; ns++) {
      if (solids[j].na == 1) {
        mass               = (silminState->fracSComp)[j][ns]*solids[j].mw;
        gibbsEnergy        = (silminState->fracSComp)[j][ns]*(solids[j].cur).g;
        enthalpy           = (silminState->fracSComp)[j][ns]*(solids[j].cur).h;
        entropy            = (silminState->fracSComp)[j][ns]*(solids[j].cur).s;
        volume             = (silminState->fracSComp)[j][ns]*(solids[j].cur).v;
        heatCapacity       = (silminState->fracSComp)[j][ns]*(solids[j].cur).cp;
        totalMass         += (silminState->fracSComp)[j][ns]*solids[j].mw;
        totalGibbsEnergy  += (silminState->fracSComp)[j][ns]*(solids[j].cur).g;
        totalEnthalpy     += (silminState->fracSComp)[j][ns]*(solids[j].cur).h;
        totalEntropy      += (silminState->fracSComp)[j][ns]*(solids[j].cur).s;
        totalVolume       += (silminState->fracSComp)[j][ns]*(solids[j].cur).v;
        totalHeatCapacity += (silminState->fracSComp)[j][ns]*(solids[j].cur).cp;
        fprintf(output, "\n%-15.15s  mass = %.2f (gm)  density = %.2f (gm/cc)", solids[j].label, mass, (volume == 0.0) ? 0.0 : mass/(10.0*volume));
        fprintf(output, "\n");
        fprintf(output, "                 %s\n", solids[j].formula); 
        fprintf(output, "                 "); 
        fprintf(output, "G = %.2f (J)  ", gibbsEnergy);
        fprintf(output, "H = %.2f (J)  ", enthalpy);
        fprintf(output, "S = %.2f (J/K)  ", entropy);
        fprintf(output, "V = %.2f (cc)  ", volume*10.0);
        fprintf(output, "Cp = %.2f (J/K)  ", heatCapacity);
        fprintf(output, "\n");
      } else {
        char *formula;
        for (i=0, mass=0.0; i<solids[j].na; i++) {
          m[i] = (silminState->fracSComp)[j+1+i][ns];
          mass += (silminState->fracSComp)[j+1+i][ns]*solids[j+1+i].mw;
        }
        (*solids[j].convert)(SECOND, THIRD, silminState->T, silminState->P, NULL, m, r, NULL, NULL, NULL, NULL, NULL);
        (*solids[j].display)(FIRST, silminState->T, silminState->P, r, &formula);
        (*solids[j].gmix) (FIRST, silminState->T, silminState->P, r, &gibbsEnergy,  NULL, NULL, NULL);
        (*solids[j].hmix) (FIRST, silminState->T, silminState->P, r, &enthalpy);
        (*solids[j].smix) (FIRST, silminState->T, silminState->P, r, &entropy,      NULL, NULL);
        (*solids[j].vmix) (FIRST, silminState->T, silminState->P, r, &volume,       NULL, NULL, NULL, NULL, NULL, NULL,  NULL, NULL, NULL);
        (*solids[j].cpmix)(FIRST, silminState->T, silminState->P, r, &heatCapacity, NULL, NULL);
        gibbsEnergy  *= (silminState->fracSComp)[j][ns]; 
        enthalpy     *= (silminState->fracSComp)[j][ns]; 
        entropy      *= (silminState->fracSComp)[j][ns];
        volume       *= (silminState->fracSComp)[j][ns];
        heatCapacity *= (silminState->fracSComp)[j][ns];
        for (i=0; i<solids[j].na; i++) {
          gibbsEnergy  += m[i]*(solids[j+1+i].cur).g;
          enthalpy     += m[i]*(solids[j+1+i].cur).h;
          entropy      += m[i]*(solids[j+1+i].cur).s;
          volume       += m[i]*(solids[j+1+i].cur).v;
          heatCapacity += m[i]*(solids[j+1+i].cur).cp;
        }
        totalMass         += mass;
        totalGibbsEnergy  += gibbsEnergy;
        totalEnthalpy     += enthalpy;
        totalEntropy      += entropy;
        totalVolume       += volume;
        totalHeatCapacity += heatCapacity;
        fprintf(output, "\n%-15.15s  mass = %.2f (gm)  density = %.2f (gm/cc)", solids[j].label, mass, (volume == 0.0) ? 0.0 : mass/(10.0*volume));
        fprintf(output, "     (analysis in mole %%)\n");
        fprintf(output, "                 %s\n", formula); 
        fprintf(output, "                 "); 
        fprintf(output, "G = %.2f (J)  ", gibbsEnergy);
        fprintf(output, "H = %.2f (J)  ", enthalpy);
        fprintf(output, "S = %.2f (J/K)  ", entropy);
        fprintf(output, "V = %.2f (cc)  ", volume*10.0);
        fprintf(output, "Cp = %.2f (J/K)  ", heatCapacity);
        fprintf(output, "\n");
        fprintf(output, "     "); 
        for (i=0; i<solids[j].na; i++) fprintf(output, " %13.13s", solids[j+1+i].label);
        fprintf(output, "\n");
        fprintf(output, "     "); 
        for (i=0; i<solids[j].na; i++) fprintf(output, " %13.2f", 100.0*(silminState->fracSComp)[j+1+i][ns]/(silminState->fracSComp)[j][ns]);
        fprintf(output, "\n");
      }
    }
  }
  
  if (silminState->fractionateLiq) {
    char *formula;
    for (i=0, mass=0.0, moles=0.0; i<nlc; i++) {
      double mw;
      for (j=0, mw = 0.0; j<nc; j++) mw += (liquid[i].liqToOx)[j]*bulkSystem[j].mw;
      m[i]   = (silminState->fracLComp)[i];
      moles += m[i];
      mass  += (silminState->fracLComp)[i]*mw;
    }
    if (mass > 0.0) {
      conLiq  (SECOND, THIRD, silminState->T, silminState->P, NULL, m, r, NULL, NULL, NULL, NULL);
      dispLiq (FIRST, silminState->T, silminState->P, r, &formula);
      gmixLiq (FIRST, silminState->T, silminState->P, r, &gibbsEnergy, NULL, NULL);
      hmixLiq (FIRST, silminState->T, silminState->P, r, &enthalpy, NULL);
      smixLiq (FIRST, silminState->T, silminState->P, r, &entropy, NULL, NULL, NULL);
      vmixLiq (FIRST, silminState->T, silminState->P, r, &volume, NULL, NULL, NULL, NULL, NULL, NULL,  NULL, NULL, NULL);
      cpmixLiq(FIRST, silminState->T, silminState->P, r, &heatCapacity, NULL, NULL);
      gibbsEnergy  *= moles;
      enthalpy     *= moles;
      entropy	   *= moles;
      volume	   *= moles;
      heatCapacity *= moles;
      for (i=0; i<nlc; i++) {
        gibbsEnergy  += m[i]*(liquid[i].cur).g;
        enthalpy     += m[i]*(liquid[i].cur).h;
        entropy      += m[i]*(liquid[i].cur).s;
        volume       += m[i]*(liquid[i].cur).v;
        heatCapacity += m[i]*(liquid[i].cur).cp;
      }
      totalMass 	+= mass;
      totalGibbsEnergy  += gibbsEnergy;
      totalEnthalpy	+= enthalpy;
      totalEntropy	+= entropy;
      totalVolume	+= volume;
      totalHeatCapacity += heatCapacity;
      fprintf(output, "\n%-15.15s  mass = %.2f (gm)  density = %.2f (gm/cc)", "liquid", mass, (volume == 0.0) ? 0.0 : mass/(10.0*volume));
      fprintf(output, "     (analysis in wt %%)\n");
      fprintf(output, " 		%s\n", formula);
      fprintf(output, " 		");
      fprintf(output, "G = %.2f (J)  ", gibbsEnergy);
      fprintf(output, "H = %.2f (J)  ", enthalpy);
      fprintf(output, "S = %.2f (J/K)  ", entropy);
      fprintf(output, "V = %.2f (cc)  ", volume*10.0);
      fprintf(output, "Cp = %.2f (J/K)  ", heatCapacity);
      fprintf(output, "\n");
    }
  }

  if (vLiq > totalVolume) fprintf(output,"\nViscosity of the System: %.2f (log 10 poise)\n", viscosity - 2.0*log10(1.0-2.0*totalVolume/(totalVolume+vLiq)));
  else fprintf(output,"\nViscosity of the System cannot be computed.\n");

  fprintf(output, "\n%-15.15s  mass = %.2f (gm)  density = %.2f (gm/cc)\n", "System", mLiq+totalMass, (totalVolume+vLiq == 0.0) ? 0.0 : 
    (mLiq+totalMass)/(10.0*(vLiq+totalVolume)));
  fprintf(output, "                 "); 
  fprintf(output, "G = %.2f (J)  ",    gLiq+totalGibbsEnergy);
  fprintf(output, "H = %.2f (J)  ",    hLiq+totalEnthalpy);
  fprintf(output, "S = %.2f (J/K)  ",  sLiq+totalEntropy);
  fprintf(output, "V = %.2f (cc)  ",   (vLiq+totalVolume)*10.0);
  fprintf(output, "Cp = %.2f (J/K)  ", cpLiq+totalHeatCapacity);
  fprintf(output, "\n");

  if (silminState->fo2Path != FO2_NONE) {
    double mO2 = -silminState->oxygen;
    for (nl=0; nl<silminState->nLiquidCoexist; nl++) for (i=0; i<nlc; i++) mO2 += (oxygen.liqToOx)[i]*(silminState->liquidComp)[nl][i];
    for (i=0; i<npc; i++) for (ns=0; ns<(silminState->nSolidCoexist)[i]; ns++) {
      if (solids[i].na == 1) mO2 += (oxygen.solToOx)[i]*(silminState->solidComp)[i][ns];
      else {
        for (j=0; j<solids[i].na; j++) mO2 += (oxygen.solToOx)[i+1+j]*(silminState->solidComp)[i+1+j][ns];
      }
    }
    fprintf(output, "\n%-15.15s  delta moles = %.6g  delta grams = %.6g\n", "Oxygen", mO2, mO2*31.9988);
    fprintf(output, "                 "); 
    fprintf(output, "G = %.2f (J)  ",    mO2*(oxygen.cur).g);
    fprintf(output, "H = %.2f (J)  ",    mO2*(oxygen.cur).h);
    fprintf(output, "S = %.2f (J/K)  ",  mO2*(oxygen.cur).s);
    fprintf(output, "V = %.2f (cc)  ",   mO2*10.0*(oxygen.cur).v);
    fprintf(output, "Cp = %.2f (J/K)  ", mO2*(oxygen.cur).cp);
    fprintf(output, "\n");
  }

  if (silminState->assimilate) {
    fprintf(output, "\nSummary of assimilant: ");
    fprintf(output, "(total mass = %.2f grams, temperature = %.2f)\n", silminState->assimMass, silminState->assimT);
    for (j=0; j<npc; j++) if (solids[j].type == PHASE)
    for (ns=0; ns<(silminState->nAssimComp)[j]; ns++) {
      if (solids[j].na == 1) {
        mass = (silminState->assimComp)[j][ns]*solids[j].mw*silminState->assimMass/silminState->dspAssimMass;
        fprintf(output, "\n%-15.15s  mass = %.2f (gm)", solids[j].label, mass);
        fprintf(output, "     (analysis in mole %%)\n");
      } else {
        for (i=0, mass=0.0; i<solids[j].na; i++) mass += (silminState->assimComp)[j+1+i][ns]*solids[j+1+i].mw;
        mass *= silminState->assimMass/silminState->dspAssimMass;
        fprintf(output, "\n%-15.15s  mass = %.2f (gm)", solids[j].label, mass);
        fprintf(output, "     (analysis in mole %%)\n");
        fprintf(output, "     "); 
        for (i=0; i<solids[j].na; i++) fprintf(output, " %13.13s", solids[j+1+i].label);
        fprintf(output, "\n");
        fprintf(output, "     "); 
        for (i=0; i<solids[j].na; i++) fprintf(output, " %13.2f", 100.0*(silminState->assimComp)[j+1+i][ns]/(silminState->assimComp)[j][ns]);
        fprintf(output, "\n");
      }
    }
  }


}