class gasState {
  double t,r,p,e,s,h,a;
  int units;
  boolean error;

  gasState(double tt,double rr, double pp, double ee, double ss, double hh, double aa, int uu) {
	t=tt;r=rr;p=pp;e=ee;s=ss;h=hh;a=aa;units=uu;error=false;
  }

  void toSI() {
    switch (units) {
      case 1:   //SI already
      break;
      case 2:   //British
      p*=6894.747;r*=16.01846;t/=1.8;e*=2325.965;
      h*=2325.965;s*=(1.8*2325.965);a*=.3048;
      units=1;
      break;
    }
  }

  void toBrit() {
    switch (units) {
      case 1:   //SI
      p/=6894.747;r/=16.01846;t*=1.8;e/=2325.965;
      h/=2325.965;s/=(1.8*2325.965);a/=.3048;
      units=2;
      break;
      case 2:   //British already
      break;
    }
  }

  int per_hs() {           //Determines p, e, and rho iteratively as functions of h and s (old airtbn routine)
      double s0,p0,t0,r0,h0,gascon,ppn,hn=0.0,affn,ffn;
      double pp[]=new double[51];
      double ff[]=new double[51];
      int j,l,ic;

      s0=6779.2;p0=101330.0;t0=273.15;r0=1.292;h0=156873;gascon=287.06;
      pp[1]=Math.exp(-(s-s0)/gascon)*Math.pow(h/h0,5.);
      pp[2]=1.e-6*pp[1];

      for(ic=1;ic<=20;ic++) {        //Bracket solution;
        for(j=1;j<=2;j++) {
            p=pp[j]*p0;
            e_ps();r_ps();
            hn=e+p/r;
            ff[j]=(hn-h)/h;
        }
        if(ff[1]*ff[2]>0.0) {pp[1]*=2.;pp[2]*=0.5;}
        else break;
      }
      if(ic>=20) {error=true;return(0);}

      for(j=1;j<=50;j++) {           //Regular Falsi
        ppn=(pp[2]*ff[1]-pp[1]*ff[2])/(ff[1]-ff[2]);
        p=ppn*p0;
        e_ps();r_ps();
        hn=e+p/r;
        ffn=(hn-h)/h;
        affn=Math.abs(ffn);
        if(affn<0.1) break;
        if(ffn*ff[1]>0.0) {pp[1]=ppn;ff[1]=ffn;}
        else {pp[2]=ppn;ff[2]=ffn;}
      }

      for(j=3;j<=50;j++) {             //Secant Method
        pp[j]=pp[j-1]-ff[j-1]*(pp[j-1]-pp[j-2])/(ff[j-1]-ff[j-2]);
        p=pp[j]*p0;
        e_ps();r_ps();
        hn=e+p/r;
        ff[j]=(hn-h)/h;
        if(Math.abs(ff[j])<5.e-4) break;
      }
      error=false;
      return(0);
}

  int r_pt() {           //Determines density iteratively as functions of p and t (old airtbn routine)
      double s0,p0,t0,r0,h0,gascon,rrn,aggn,ggn,temp;
      double rr[]=new double[51];
      double gg[]=new double[51];
      int j,l,ic;

      s0=6779.2;p0=101330.0;t0=273.15;r0=1.292;h0=156873;gascon=287.06;
      rr[1]=p*t0/t/p0;
      rr[2]=0.9*rr[1];
      temp=t;

      for(ic=1;ic<=10;ic++) {        //Bracket solution;
        for(j=1;j<=2;j++) {
            r=rr[j]*r0;
            t_pr();
            gg[j]=(t-temp)/temp;
        }
        if(gg[1]*gg[2]>0.0) {rr[1]*=2.;rr[2]*=0.5;}
        else break;
      }
      if(ic>=10) {error=true;return(0);}

      for(j=1;j<=50;j++) {           //Regular Falsi
        rrn=(rr[2]*gg[1]-rr[1]*gg[2])/(gg[1]-gg[2]);
        r=rrn*r0;
        t_pr();
        ggn=(t-temp)/temp;
        aggn=Math.abs(ggn);
        if(aggn<0.1) break;
        if(ggn*gg[1]>0.0) {rr[1]=rrn;gg[1]=ggn;}
        else {rr[2]=rrn;gg[2]=ggn;}
      }

      for(j=3;j<=50;j++) {             //Secant Method
        rr[j]=rr[j-1]-gg[j-1]*(rr[j-1]-rr[j-2])/(gg[j-1]-gg[j-2]);
        r=rr[j]*r0;
        t_pr();
        gg[j]=(t-temp)/temp;
        if(Math.abs(gg[j])<5.e-4) break;
      }
      t=temp;error=false;
      return(0);
  }

  int r_ps() {                     //Determines density as a function or pressure and entropy  (old tgas5.f)
      double p0 = 101330.0,r0 = 1.292,s0 = 6779.2,gascon = 287.6;
      double x,snon,z,delts,rratio,y,zm;

      x = Math.log(p/p0)/Math.log(10);
      snon = Math.log(s/gascon)/Math.log(10.);
      z = x - snon;

      if (snon < 1.23) {
        delts = s - s0;
        rratio = (Math.log(p/p0))/(Math.log(10.)) / 1.4-delts / (3.5*gascon);
        r = r0*Math.exp(rratio);
        return(0);
      }
      if (snon < 1.42) {              // line 10
        y=-17.2119+54.9354*snon;
        y+=(-1.99776+3.17884*snon)*z;
        y+=(-46.9831-0.866580*z+12.1069*snon)*snon*snon;
        y+=(0.158567-0.103055*snon-0.00152322*z)*z*z;
        r=Math.pow(10.,y)*r0;
        return(0);
      }
      if (snon < 1.592 ) {  // line 20
        y=280.393-595.834*snon+(-17.5934+24.9706*snon)*z+(421.767-8.85461*z-99.4515*snon)*snon*snon+(0.36189-0.255458*snon-0.00296892*z)*z*z;
        y/=1.0+Math.exp(-15.0*(x+54.179-86.947*snon+33.583*snon*snon));
        y+=-278.074+611.791*snon+(13.7528-19.2394*snon)*z+(-442.909+7.10425*z+105.869*snon)*snon*snon+(-0.197269+0.149708*snon+0.00119153*z)*z*z;
        r=Math.pow(10,y)*r0;
        return(0);
      }
      if (snon<1.70) { // line 30
        zm = 7.5269*snon-14.9366;
        if (z>=zm) { // line 40
           y=-50.1859+63.1564*snon+(2.39925-1.47883*snon)*z-19.8852*snon*snon+(0.466791-0.306926*snon)*z*z;
           y/=1.0+Math.exp(-2.0*(z+11.16*snon-18.203));
           y+=39.7149-48.0988*snon+(-1.17821+1.27522*snon)*z+14.5552*snon*snon+(-0.250279+0.158399*snon)*z*z;
         }
         else {
           y=-68.5292+82.3834*snon+(-1.24942+0.831615*snon)*z-24.524*snon*snon+(-0.113019+0.0746719*snon)*z*z;
           y/=1.0+Math.exp(-2.0*(z+33.976*snon-49.659));
           y+=110.732-133.968*snon+(0.37583+0.277887*snon)*z+40.3018*snon*snon+(0.118506-0.0698812*snon)*z*z;
         }
       r=Math.pow(10,y)*r0;
       return(0);
     }
     if (snon<1.80) {  // line 50
       zm = 7.5269*snon-14.9366;
       if(z<zm) {
         y=-72.8463+94.2375*snon+(4.04979-2.45467*snon)*z-30.2865*snon*snon+(0.12872-0.0863721*snon)*z*z;
         y/=1.0+Math.exp(-2.0*(z+13.4576*snon-18.36));
         y+=62.2639-78.884*snon+(-2.19123+1.84533*snon)*z+24.8413*snon*snon+(-0.0876105+0.0555047*snon)*z*z;
       }
       else {
         y=60.5677-70.8307*snon+(-2.95622+1.70504*snon)*z+20.7549*snon*snon+(-0.270962+0.180057*snon)*z*z;
         y/=1.0+Math.exp(-2.0*(z+16.822*snon-29.139));
         y+=-29.9578+38.9998*snon+(2.90256-1.25088*snon)*z-12.6641*snon*snon+(0.0613867-0.0517059*snon)*z*z;
       }
       r=Math.pow(10,y)*r0;
       return(0);
     }
     if(snon<1.90) { // line 70
       y=-639.514+697.458*snon+(-100.154+106.701*snon)*z+(-190.745-28.6323*z)*snon*snon+(-0.481471+0.200348*snon-0.00643371*z)*z*z;
       y/=1.0+Math.exp(-2.0*(z+35.275*snon-58.624));
       y+=623.124-677.571*snon+(91.2811-95.548*snon)*z+(184.603+25.4274*z)*snon*snon+(0.697635-0.327916*snon+0.00490838*z)*z*z;
       r=Math.pow(10,y)*r0;
       return(0);
     }
     if(snon<2.00) {      // line 80
       y=-113.217+124.304*snon+(4.41505-2.4853*snon)*z-34.2370*snon*snon+(-0.0784297+0.0121459*snon)*z*z;
       y/=1.0+Math.exp(-2.00*(z+20.884*snon-36.54));
       y+=14.0088-15.8855*snon+(0.171245+0.432352*snon)*z+4.42836*snon*snon+(-0.00954417+0.00892335*snon)*z*z;
       r=Math.pow(10,y)*r0;
       return(0);
     }
     if(snon<=2.10) {   // line 90
       y=(65.2920-62.8154*snon)+(2.10906-1.11759*snon)*z+15.0026*snon*snon+(0.271716-0.14431*snon)*z*z;
       y/=1.0+Math.exp(-5.00*(z+28.284*snon-53.185));
       y+=-36.2767+38.634*snon+(1.48507-0.341824*snon)*z-10.4718*snon*snon+(-0.0325437+0.0122032*snon)*z*z;
       r=Math.pow(10,y)*r0;
       return(0);
     }
     if(snon>2.10) {  // line 100
       y=-127.867+85.8453*snon+(-22.8808+11.3034*snon)*z-11.234*snon*snon+(-1.56005+0.779397*snon)*z*z;
       y/=1.0+Math.exp(-2.00*(z+15.048*snon-26.307));
       y+=186.938-145.261*snon+(22.3883-10.4974*snon)*z+26.0284*snon*snon+(1.85949-0.898218*snon)*z*z;
       r=Math.pow(10,y)*r0;
       return(0);
     }
     return(0);
  }

  int a_ps() {	// speed of sound as a function of pressure and entropy (old tgas7.f)

	double p0 = 101330,	s0 = 6779.2, r0 = 1.292, a0 = 331.3613, gascon = 287.06;
	double x,snon,y,z,zm,delts,asqlog;

	x=(Math.log(p/p0))/(Math.log(10));
	snon=(Math.log(s/gascon))/(Math.log(10));
	z = x - snon;

	if(snon < 1.23) {
		delts = s-s0;
		asqlog = Math.log(1.4*p0/r0)+(Math.log(p/p0)+delts/gascon)/3.5;
		a = Math.exp(asqlog/2.0);
		return(0);}
	else if(snon <= 1.4) y = (-.138377-8.84138*snon) + ((2.6105-3.16535*snon)*z) + ((11.0866+.988389*z-3.25761*snon)*snon*snon) + ((-.100224+.0662193*snon+.00082061*z)*z*z);
	else if(snon < 1.595) y = (131.057-288.847*snon) + ((-5.04887+7.73862*snon)*z) + ((210.147-2.88963*z-50.396*snon)*snon*snon) + ((.0548031-.0439459*snon-.000210202*z)*z*z) + ( (-133.465+284.739*snon) + ((7.57389-10.7749*snon)*z) + ((-202.362+3.8313*z+47.9075*snon)*snon*snon) + ((-.153453+.108531*snon+.00097931*z)*z*z))/(1.0+(Math.exp(-15.0*(x+54.179-86.947*snon+33.583*snon*snon))));
	else if(snon < 1.693) {
		zm = -9.842*snon+14.19;
		if(z <= zm) y = (-61.3548+78.0742*snon) + ((2.08524-1.21609*snon)*z) + (-24.3686*snon*snon) + ((.0877563-.0546311*snon)*z*z) + ( (20.7952-27.1591*snon) + ((-.743673+.408312*snon)*z) + (8.68124*snon*snon) + ((-.093592+.0532328*snon)*z*z))/(1.0+(Math.exp(-2.0*(z+38.785*snon-57.157))));
		else y = (3.37056-4.87016*snon) + ((-.385754+.287192*snon)*z) + (2.02041*snon*snon) + ((-.00463144+.00830832*snon)*z*z);}
	else if(snon <1.80) {
		zm = -1.917*snon+0.092;
		if( z <= zm) y = (-80.4927+76.3739*snon) + ((-9.381+5.72104*snon)*z) + ( -16.3435*snon*snon) + ((-.748578+.450043*snon)*z*z) + ((83.4054-85.8837*snon) + ((4.84197-3.11188*snon)*z) + (21.1196*snon*snon) + ((.233945-.159099*snon)*z*z))/(1.0+(Math.exp(-2.0*(z+7.874*snon-7.569))));
		else y = (-4.73308+4.69363*snon) + ((.0943798+.00354953*snon)*z) + (-.798293*snon*snon) + ((.0202561-.00873036*snon)*z*z);}
	else if(snon < 1.9) y = (-660.574+738.042*snon) + ((-87.7589+97.894*snon)*z) + ((-206.156-27.3753*z)*snon*snon) + (-.0214028*z*z) + ((665.014-743.416*snon) + ((88.7679-99.0508*snon)*z) + ((208.117+27.7374*z)*snon*snon) + (.00281148*z*z))/(1.0+(Math.exp(-2.0*(z+35.275*snon-58.624))));
	else if(snon < 2.0) y = (-5.93554-7.79929*snon) + ((-7.23618+3.31162*snon)*z) + (5.06381*snon*snon) + ((-.5377735+.246865*snon)*z*z) + ((39.526-29.0994*snon) + ((6.24136-2.70007*snon)*z) + (5.42786*snon*snon) + ((.595235-.267771*snon)*z*z))/(1.0+(Math.exp(-2.0*(z+13.6751*snon-20.1676))));
	else if(snon < 2.1) y = (70.2453-73.4732*snon) + (-.586844*z) + ((19.5548+.181521*z)*snon*snon) + (.00736786*z*z) + ((-16.961+18.0169*snon) + ( -.365779*z) + ((-4.75063+.107583*z)*snon*snon) + (.00201803*z*z))/(1.0+(Math.exp(-5.0*(z+28.284*snon-53.185))));
	else y = (-112.793+135.296*snon) + (2.74295*z) + ((-40.8171-1.29257*z)*snon*snon) + (-.225652*z*z) + ((257.057-283.199*snon) + ( -5.84656*z) + ((79.1574+2.17312*z)*snon*snon) + (.335658*z*z))/(1.0+(Math.exp(-2.0*(z+15.048*snon-26.307))));

	a = Math.pow(10,y)*a0;
	return(0);
  }

  int e_ps() { // energy as a function of pressure and entropy (old tgas6.f)

	double p0=101330,e0 = 78408.4,s0 = 6779.2,gascon = 287.06;
	double z,snon,x,delts,eratio,y,zm;

	x=(Math.log(p/p0))/(Math.log(10));
	snon=(Math.log(s/gascon))/(Math.log(10));
	z = x - snon;

	if (snon < 1.23) {
		delts = s - s0;
		eratio = ((Math.log(p/p0))/(Math.log(10))+delts/gascon)/3.5;
		e = 2.5*e0*Math.exp(eratio);
		return(0);}
	else if(snon < 1.4) y = (48.9511-115.989*snon) + ((-.443026+1.49819*snon)*z) + ((87.0709-.702417*z-20.4916*snon)*snon*snon) + ((-.0509747+.0294929*snon+.000524685*z)*z*z);
	else if(snon < 1.592) y = ((130.563-287.254*snon)) + ((-7.71873+1.17605*Math.pow(10,1)*snon)*z) + ((206.757-4.31838*z-48.3472*snon)*snon*snon) + ((.164299-.121768*snon-.00126118*z)*z*z) + ((-84.598+174.701*snon) + ((9.52234-13.2819*snon)*z) + ((-120.234+4.6271*z+27.5818*snon)*snon*snon) + ((-.227725+.160365*snon+.0014933*z)*z*z)) /(1.0 + (Math.exp(-15.0*(x+54.179-86.947*snon+33.583*snon*snon))));
	else if(snon < 1.7)  {
		zm = -1.917*snon+0.092;
		if(z <= zm) y = (-134.875+192.654*snon) + ((10.4478-6.60355*snon)*z) + (-67.2013*snon*snon) + ((.414546-.275698*snon)*z*z) + ((107.768-161.662*snon) + ((-10.1257+6.4962*snon)*z) + (59.3194*snon*snon) + ((-.463204+.304448*snon)*z*z)) /(1.0 + (Math.exp(-2.0*(z-50.04*snon-74.698))));
		else y = (-33.9933+40.2792*snon) + ((.762742-.379044*snon)*z) + (-10.9895*snon*snon) + ((-.0713412+.0422359*snon)*z*z) + ((33.2182-40.3847*snon) + ((.187601-.107599*snon)*z) + (12.2657*snon*snon) + ((-.165537+.10256*snon)*z*z)) /(1.0 + (Math.exp(-2.0*(z+14.062*snon-24.046))));}
	else if(snon < 1.8) {
		zm = -1.917*snon+0.092;
		if(z <= zm) y = (-16.6572+15.1243*snon) + ((-.472897+.295848*snon)*z) + (-2.31418*snon*snon) + ((-.0361191+.017288*snon)*z*z) + ((-3.27168+9.45383*snon) + ((2.04-1.1342*snon)*z) + (-4.28765*snon*snon) + ((.15021-.0827388*snon)*z*z)) /(1.0 + (Math.exp(-2.0*(z+16.995*snon-23.354))));
		else y = (-19.4882+23.1472*snon) + ((1.00783-.518103*snon)*z) + (-5.92829*snon*snon) + ((.0992242-.0567928*snon)*z*z) + ((-2.09813+1.22781*snon) + ((-.627531+.371651*snon)*z) + (.000984125*snon*snon) + ((-.124024+.0765034*snon)*z*z)) /(1.0 + (Math.exp(-2.0*(z+18.152*snon-31.096))));}
	else if(snon < 1.9) y = (-63.7715+54.6005*snon) + ((-8.90073+5.02063*snon)*z) + (-9.67094*snon*snon) + ((-.713493+.397962*snon)*z*z) + ((50.5131-39.7363*snon) + ((9.27248-5.15576*snon)*z) + (6.45632*snon*snon) + ((.72563-.403054*snon)*z*z)) /(1.0 + (Math.exp(-2.0*(z+35.275*snon-58.624))));
	else if(snon < 2.0) y = (-17.5434+20.0307*snon) + ((.620484-.276756*snon)*z) + (-4.76084*snon*snon) + ((.0114921-.00580043*snon)*z*z) + ((11.0485-12.61*snon) + ((-.659043+.366266*snon)*z) + (3.59417*snon*snon) + ((-.00604061+.00678382*snon)*z*z)) /(1.0 + (Math.exp(-2.0*(z+20.884*snon-36.54))));
	else if(snon < 2.1) y = (-11.8487+12.6558*snon) + ((.000609885+.0590396*snon)*z) + (-2.45275*snon*snon) + ((-.00670888+.00521033*snon)*z*z) + ((3.39157-3.835*snon) + ((-.214031+.0990508*snon)*z) + (1.07008*snon*snon) + ((-.0787516+.0363642*snon)*z*z)) /(1.0 + (Math.exp(-5.0*(z+28.284*snon-53.185))));
	else y = (-24.0656+3.35589*snon) + ((-11.4521+5.62572*snon)*z) + (4.93842*snon*snon) + ((-.755608+.369774*snon)*z*z) + ((22.3905-1.90842*snon) + ((10.5618-5.17247*snon)*z) + (-4.3714*snon*snon) + ((.627818-.312015*snon)*z*z)) /(1.0 + (Math.exp(-2.0*(z+15.048*snon-26.307))));

	e = Math.pow(10,y)*e0;
	return(0);

  }

int h_pr() { // enthalpy as a function of pressure and density (old tgas4.f)

	double r0 = 1.292, p0 = 101330;
	double y,x,z1,hhigh=0.0,hlow,yhigh=0.0,ylow=0.0,ym=0.0,gamm=0.0,rho;
	int jflag=0,iflag=0;

	rho=r;
	y = Math.log(rho/r0)/Math.log(10);
	x = Math.log(p/p0)/Math.log(10);

    if (Math.abs(y+4.5) < .045 ) {
		iflag=0;jflag=-1;ym=y;y=-4.5+.045;yhigh=y;rho = Math.pow(10,y)*r0;
	}
	else if (Math.abs(y+0.5) < .005 ) {
		iflag=1;jflag=-1;ym=y;y=-0.5+.005;yhigh=y;rho = Math.pow(10,y)*r0;
	}
	else {iflag=-1;jflag=-1;}

    for (;;) {
	    z1=x-y;
	    if (y> -0.5) {
		    if (z1 <= 0.1) gamm=1.4017;
		    else if (z1 <= 1.05) gamm=(-96.7488+.205296*y)+((269.927-1.92887*y)*z1)+((.378392-.324965*z1-.00361036*y)*y*y)+((-246.711+1.54416*y+74.8760*z1)*z1*z1)+((98.1502-.205488*y)+((-269.913+1.93052*y)*z1)+((-.378527+.324832*z1+.00366182*y)*y*y)+((246.63-1.54646*y-74.898*z1)*z1*z1))/(1.0-(Math.exp(-26.59865+1.564631*y+23.12926*z1-1.360543*y*z1)));
		    else if (z1 <= 1.6) gamm=(-.267593-.187457*y)+((5.07693+.272286*y)*z1)+((.0104541-.0142211*z1+.000638962*y)*y*y)+((-5.0852-.0781935*y+1.58711*z1)*z1*z1)+((2.87969+.39009*y)+((-8.06179-.55125*y)*z1)+((-.0101903+.0135906*z1-.000897772*y)*y*y)+((7.29592+.183861*y-2.15153*z1)*z1*z1))/(1.0+(Math.exp(182.8573-34.28596*y-151.786*z1+29.7621*y*z1)));
		    else if (z1<=2.3) gamm=(.921537-.23967*y)+((1.30714+.34299*y)*z1)+((-.0218847+.0136691*z1-.000490274*y)*y*y)+((-1.20916-.110206*y+.308792*z1)*z1*z1)+((-6.77089-.0690476*y)+((8.18168-.0952708*y)*z1)+((.0298487-.0178706*z1+.000628419*y)*y*y)+((-3.07662+.0660408*y+.33859*z1)*z1*z1))/(1.0+(Math.exp(159.16669+39.76192*y-79.66199*z1-16.667*y*z1)));
		    else error=true;}

	    else if (y > -4.5) {
	        if (z1 <= 0.1) gamm=1.399;
		    else if (z1 <= 0.95) gamm=(-133.083-9.98707*y)+((394.734+23.581*y)*z1)+((1.43957-1.43175*z1+.0000177068*y)*y*y)+((-384.712-13.6367*y+124.325*z1)*z1*z1)+((134.486+9.99122*y)+((-394.719-23.5853*y)*z1)+((-1.43799+1.43039*z1+.000144367*y)*y*y)+((384.616+13.6318*y-124.348*z1)*z1*z1))/(1.0+(Math.exp(-21.41444+1.381584*y+20.39473*z1-1.315789*y*z1)));
		    else if (z1 <=1.5) gamm=(-7.36684-1.13247*y)+((24.7879+1.99625*y)*z1)+((-.049163+.0416673*z1-.000658149*y)*y*y)+((-23.299-.859418*y+7.19016*z1)*z1*z1)+((-2.42647+.557912*y)+((-2.03055-1.22031*y)*z1)+((.0374866-.0339278*z1+.000521042*y)*y*y)+((7.75414+.608488*y-3.68326*z1)*z1*z1))/(1.0+(Math.exp(80.77385-12.73807*y-65.47623*z1+11.90475*y*z1)));
		    else if (z1 <= 2.00) gamm=(.43152-.283857*y)+((2.27791+.399159*y)*z1)+((-.0129444+.00878724*z1-.000160583*y)*y*y)+((-1.84314-.128136*y+.445362*z1)*z1*z1)+((-10.3883-.358718*y)+((13.5068+.187268*y)*z1)+((-.00428184-.000952016*z1-.0000410506*y)*y*y)+((-5.63894-.00145626*y+.739915*z1)*z1*z1))/(1.0+(Math.exp(294.9221+13.6866*y-155.9335*z1-3.787766*y*z1)));
		    else if (z1 <= 2.5) gamm=(-3.77766-.553738*y)+((6.60834+.487181*y)*z1)+((-.0211045+.00967277*z1-.00021942*y)*y*y)+((-2.94754-.102365*y+.43962*z1)*z1*z1)+((40.5813+3.25692*y)+((-47.9583-2.5366*y)*z1)+((.0906436-.0347578*z1+.00100077*y)*y*y)+((18.904+.494114*y-2.48554*z1)*z1*z1))/(1.0+(Math.exp(534.718+74.95657*y-221.9822*z1-30.17229*y*z1)));
		    else error=true;}

 		else if ( z1 > 0.1) {
		    if (z1 <= 0.85) gamm=(253.908+101.491*y)+((-387.199-154.304*y)*z1)+((7.28532-8.04378*z1-.00182577*y)*y*y)+((98.6233+46.3763*y+21.8994*z1)*z1*z1)+((-252.423-101.445*y)+((387.21+154.298*y)*z1)+((-7.2773+8.04277*z1+.00228399*y)*y*y)+((-98.7576-46.3883*y-21.9438*z1)*z1*z1))/(1.0-(Math.exp(-11.0+2.0*y+11.0*z1-2.0*y*z1)));
		    else if (z1 <= 1.30) gamm=(-10.5745-1.9369*y)+((30.7202+3.35578*y)*z1)+((-.0779965+.066879*z1-.000986882*y)*y*y)+((-26.0637-1.42391*y+7.23223*z1)*z1*z1)+((-18.6342+.0241997*y)+((32.0880-.746914*y)*z1)+((.0375161-.0410125*z1+.000574637*y)*y*y)+((-16.9985+.539041*y+2.56253*z1)*z1*z1))/(1.0+(Math.exp(276.8567+21.52383*y-216.4837*z1-13.94837*y*z1)));
		    else if (z1 <= 1.95) gamm=(.617584-.24069*y)+((1.95904+.341644*y)*z1)+((-.0101073+.00677631*z1-.000115922*y)*y*y)+((-1.68951-.110932*y+.426058*z1)*z1*z1)+((-13.4222-.543713*y)+((18.1528+.395928*y)*z1)+((-.00741105+.00167768*z1-.00000332714*y)*y*y)+((-7.97425-.0580593*y+1.12448*z1)*z1*z1))/(1.0+(Math.exp(86.77803-8.370349*y-40.74084*z1+7.407405*y*z1)));
		    else if (z1 <= 2.6) gamm=(-8.32595-.350219*y)+((13.6455+.35935*y)*z1)+((-.00370109+.00330836*z1+.000110018*y)*y*y)+((-6.49007-.0838594*y+1.02443*z1)*z1*z1)+((-30.8441-1.4951*y)+((30.0585+.91965*y)*z1)+((-.0360024+.0102522*z1-.00046876*y)*y*y)+((-9.33522-.135228*y+.892634*z1)*z1*z1))/(1.0+(Math.exp(88.00047-16.79356*y-33.33353*z1+8.465574*y*z1)));
		    else error=true; }

		else gamm=1.3986;

        /* Out of Range */
		if(error) return(0);

		h=gamm/(gamm-1.0)*p/rho;

		/* Immediate Exit */
		if(iflag==-1 || jflag==1) return(0);

		/* Exit with Calculation */
		if(jflag==0) {hlow=h;h=hlow+(hhigh-hlow)/(yhigh-ylow)*(ym-ylow);return(0);}

		/* Iterate */
    	y=-4.5-0.045;
	    if(iflag==1) y=-0.5-0.005;
		hhigh=h;ylow=y;rho=Math.pow(10,y)*r0;jflag=0;}

    }

int s_re() {  // entropy as a function of energy and density (old tgas2.f)

	double e0 = 78408.4,r0 = 1.292,	gascon = 287.06;
	double eratio,rratio,y,z,slow,shigh=0.0,ylow=0.0,yhigh=0.0,ym=0.0,rsave=0.0,delts,deltz;
	int iflag=0,jflag=0;

	eratio = e/e0;rratio = r/r0;
	y = Math.log(rratio)/Math.log(10);
	z = Math.log(eratio)/Math.log(10);


	if(Math.abs(y+4.5)<.045) {iflag=0;jflag=-1;rsave=r;ym=y;y=-4.5+.045;yhigh=y;r=Math.pow(10,y)*r0;}
	else if(Math.abs(y+0.5)<.005){iflag=1;jflag=-1;rsave=r;ym=y;y=-0.5+0.005;yhigh=y;r=Math.pow(10,y)*r0;}
	else iflag=-1;

    for(;;) {
	  	if (z<=0.65) {
			deltz = z-0.4;
			delts=(2.5*deltz-y)*gascon*2.302585;
			s = 6779.2004+delts;}
		else if(y>-4.5) {
			if(y>-0.5) {
			    if(z>3.0) error=true;
			    s = gascon*( (20.1858-3.13458*y) + ((10.3619+1.87767*y)*z) + ((-.172922+.112174*z+.012826*y)*y*y) + ((-5.43557-.871048*y+2.01789*z)*z*z) );}
			else {
			    if(z>3.4) error=true;
		    	s = gascon*((10.836-4.55524*y) + ((29.6473+3.90851*y)*z) + ((-.00205732+.0365982*z+.00523821*y)*y*y) + ((-16.7001-1.44623*y+3.98307*z)*z*z));}}
        else {
            if(z>3.69) error=true;
		    s = gascon*((-.991081-5.00277*y) + ((54.6521+5.10144*y)*z) + ((.0176206+.0212002*z+.00176358*y)*y*y) + ((-29.7001-1.84915*y+5.87892*z)*z*z));}

	    if(iflag==-1 || jflag==1) return(0);
	    if(jflag==0) {
			slow=s;
			s=slow+(shigh-slow)/(yhigh-ylow)*(ym-ylow);
			r=rsave;
			return(0);}
	    shigh=s;
	    y=-4.5-.045;
	    if(iflag==1) y = -0.5+.005;
	    ylow = y;r = Math.pow(10,y)*r0;jflag=0;}

    }

int pa_er(int kflag) {  // pressure and speed of sound as a function of energy and density (old tgas.f)
                        // kflag=0 for just pressure, kflag=1 for pressure and sound speed
	int iflag=0,jflag=0;
	double e0 = 78408.4, r0 = 1.292;
	double rratio,eratio,y,z,yn,yhigh=0.,ylow=0.,rsave=0.,gamm,gammr=0.,gamme=0.,phigh=0.,plow,ahigh=0.,alow,ym=0.,deno;
	double[] g=new double[10];
	double[] gr=new double[10];
	double[] ge=new double[10];

	rratio = r/r0;
	eratio = e/e0;
	y = Math.log(rratio)/Math.log(10);
	z = Math.log(eratio)/Math.log(10);

    if(Math.abs(y+4.5)<.025) {iflag=0;jflag=-1;rsave=r;ym=y;y=-4.5+.025;yhigh=y;r=Math.pow(10,y)*r0;}
	else if(Math.abs(y+0.5)<.005){iflag=1;jflag=-1;rsave=r;ym=y;y=-0.5+.005;yhigh=y;r=Math.pow(10,y)*r0;}
    else iflag=-1;


	for(;;) {
		if(y>-0.5) {
			if(z<=0.65) {gamm=1.396;if(kflag!=0) {gammr=0.0;gamme=0.0;}}
			else if(z<=1.7) {
				g[1] = 1.37062+.0129673*y;
				g[2] = (.111418-.0326912*y)*z;
				g[3] = (.00106869-.00200286*z+.000238305*y)*y*y;
				g[4] = (-.106133+.0190251*y+.0030221*z)*z*z;
				gamm=g[1]+g[2]+g[3]+g[4];
				if(kflag!=0) {
    				gr[1] = .0129673;
	    			gr[2] = -.0326912*z;
		    		gr[3] = (.00213738-.00400572*z+.000714915*y)*y;
			    	gr[4] = .0190251*z*z;
				    ge[2] = g[2]/z;
		    		ge[3] = -.00200286*y*y;
		    		ge[4] = (-.212266+.0380502*y+.0090663*z)*z;
			    	gammr = gr[1]+gr[2]+gr[3]+gr[4];
			    	gamme=ge[2]+ge[3]+ge[4];}}
			else if(z<=2.35) {
				g[1] = .0343846-.233584*y;
				g[2] = (2.85574+.259787*y)*z;
				g[3] = (-.01089927+.00423659*z+.000385712*y)*y*y;
				g[4] = (-1.94785-.0673865*y+.408518*z)*z*z;
				g[5] = -4.20569+.133139*y;
				g[6] = (4.51236-.166341*y)*z;
				g[7] = (.00167787-.00110022*z+.000306676*y)*y*y;
				g[8] = (-1.35516+.0491716*y+.0752509*z)*z*z;
				g[9] = (175.7042-2.163278*y-88.33702*z+1.897543*y*z);
				if(kflag!=0) {
				    gr[1] =-.0233584;
			    	gr[2] = .259787*z;
			    	gr[3] = (-.02179854+.00847318*z+.001157136*y)*y;
			    	gr[4] = -.0673865*z*z;
			    	gr[5] = .133139;
			    	gr[6] = -.166341*z;
			    	gr[7] = (.00335574-.00220044*z+.000920028*y)*y;
			    	gr[8] = .0491716*z*z;
			    	gr[9] = -2.163278+1.897543*z;
			    	ge[2] = g[2]/z;
			    	ge[3] = .00423659*y*y;
			    	ge[4] = (-3.8957-.134773*y+1.225554*z)*z;
			    	ge[6] = g[6]/z;
			    	ge[7] = -.00110022*y*y;
			    	ge[8] = (-2.71032+.0983432*y+.2257527*z)*z;
			    	ge[9] = -88.33702+1.897543*y;}
				g[9] = Math.exp(g[9]);	// line 210
				gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
				if (kflag!=0) {
   					gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
    				gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}
			else {
		        if(z>2.9) error=true;
				g[1] = -1.70633-.148403*y;
				g[2] = (4.23104+.13729*y)*z;
				g[3] = (-.00910934+.00385707*z+.000269026*y)*y*y;
				g[4] = (-1.97292-.028183*y+.295882*z)*z*z;
				g[5] = 34.158-18.9972*y;
				g[6] = (-40.858+13.0321*y)*z;
				g[7] = (-.801272+.275121*z-.000177969*y)*y*y;
				g[8] = (16.0826-2.23386*y-2.08853*z)*z*z;
				g[9] = (256.1323+173.7089*y-90.5889*z-58.38803*y*z);
				if (g[9]>30.0) g[9] = 30.0;
				if (g[9]<-30.0) g[9] = -30.0;
				if(kflag!=0) {
				    gr[1] = -.148403;
				    gr[2] = .13729*z;
				    gr[3] = (-.01821868+.00771414*z+.000807078*y)*y;
				    gr[4] = -.028183*z*z;
				    gr[5] = -18.9972;
	    		    gr[6] = 13.0321*z;
				    gr[7] = (-1.602544+.550242*z-.000533907*y)*y;
				    gr[8] = -2.23386*z*z;
				    gr[9] = 173.7089-58.38803*z;
				    ge[2] = g[2]/z;
				    ge[3] = .00385707*y*y;
				    ge[4] = (-3.94584-.056366*y+.887646*z)*z;
				    ge[6] = g[6]/z;
				    ge[7] = .275121*y*y;
				    ge[8] = (32.1652-4.46772*y-6.26559*z)*z;
				    ge[9] = -90.5889-58.38803*y;}
				g[9] = Math.exp(g[9]);	// line 210
				gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
				if (kflag!=0) {
    				gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
	    			gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}}
		else if(y>-4.5) {
			if(z<=0.65) {gamm=1.396;if(kflag!=0) {gammr=0.0;gamme=0.0;}}
			else if(z<=1.5) {
				g[1] = 1.39123-.00408321*y;
				g[2] = (.0142545+.0141769*y)*z;
				g[3] = (.000257225+.000652912*z+.0000846912*y)*y*y;
				g[4] = (.062555-.00783637*y-.097872*z)*z*z;
				g[5] = 5.80955-.182302*y;
				g[6] = (-9.62396+.179619*y)*z;
				g[7] = (-.0230518+.011872*z-.000335499*y)*y*y;
				g[8] = (5.27047-.0365507*y-.919897*z)*z*z;
				g[9] = (-10.0*z+14.2);
				if (kflag!=0) {
    				gr[1] = -.00408321;
	    			gr[2] = .0141769*z;
		    		gr[3] = (.00051445+.001305824*z+.0002540736*y)*y;
			    	gr[4] = -.00783637*z*z;
				    gr[5] = -.182302;
    				gr[6] = .179619*z;
	    			gr[7] = (-.0461036+.023744*z-.001006497*y)*y;
		    		gr[8] = -.0365507*z*z;
			    	gr[9] = 0.0;
			    	ge[2] = g[2]/z;
			    	ge[3] = .000652912*y*y;
			    	ge[4] = (.12511-.01567274*y-.293616*z)*z;
			    	ge[6] = g[6]/z;
			    	ge[7] = .011872*y*y;
			    	ge[8] = (10.54094-.0731014*y-2.759691*z)*z;
			    	ge[9] = -10.0;}
			    g[9] = Math.exp(g[9]);	// line 210
			    gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
			    if (kflag!=0) {
   				    gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
    		 	    gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}
			else if(z<=2.22) {
				g[1] = -1.20784-.257909*y;
				g[2] = (5.02307+.287201*y)*z;
				g[3] = (-.00995577+.00523524*z-.000145574*y)*y*y;
				g[4] = (-3.20619-.0750405*y+.651564*z)*z*z;
				g[5] = -6.62841+.0277112*y;
				g[6] = (7.30762-.076823*y)*z;
				g[7] = (.00719421-.00362463*z+.000162777*y)*y*y;
				g[8] = (-2.33161+.0304767*y+.166856*z)*z*z;
				g[9] = (125.5324+2.015335*y-63.90747*z-.6515225*y*z);
				if (kflag!=0) {
					gr[1] = -.257909;
    				gr[2] = .287201*z;
	    			gr[3] = (-.01991154+.01047048*z-.000436722*y)*y;
		    		gr[4] = -.0750405*z*z;
			    	gr[5] = .0277112;
					gr[6] = -.076823*z;
    				gr[7] = (.01438842-.00724926*z+.000488331*y)*y;
	    			gr[8] = .0304767*z*z;
		    		gr[9] = 2.015335-.06515225*z;
			    	ge[2] = g[2]/z;
				    ge[3] = .00523524*y*y;
					ge[4] = (-6.41238-.150081*y+1.954692*z)*z;
    				ge[6] = g[6]/z;
	    			ge[7] = -.00362463*y*y;
		    		ge[8] = (-4.66322+.0609534*y+.500568*z)*z;
			    	ge[9] = -63.90747-.6515225*y;}
			    g[9] = Math.exp(g[9]);	// line 210
			    gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
			    if (kflag!=0) {
				    gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
    		 	    gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}
            else if(z<=2.95) {
    			g[1] = -2.2646-.0782263*y;
	    		g[2] = (4.90497+.0718096*y)*z;
	    		g[3] = (-.00306443+.00174209*z+.0000284214*y)*y*y;
	    		g[4] = (-2.2475-.0131641*y+.333658*z)*z*z;
	    		g[5] = -14.7904-.176627*y;
	    		g[6] = (13.5036+.087728*y)*z;
    			g[7] = (-.00213327+.000715487*z+.0000730928*y)*y*y;
    			g[8] = (-3.95372-.00896151*y+.363229*z)*z*z;
    			g[9] = (178.8542+6.317894*y-67.56741*z-2.46006*y*z);
		    	if (kflag!=0) {
	                gr[1] = -.0782263;
        			gr[2] = .0718096*z;
        			gr[3] = (-.00612886+.00348418*z+.0000852642*y)*y;
        			gr[4] = -.0131641*z*z;
        			gr[5] = -.176627;
        			gr[6] = .087728*z;
        			gr[7] = (-.00426654+.001430974*z+.0002192784*y)*y;
        			gr[8] = -.00896151*z*z;
        			gr[9] = 6.317894-2.46006*z;
        			ge[2] = g[2]/z;
        			ge[3] = .00174209*y*y;
           			ge[4] = (-4.495-.0263282*y+1.000974*z)*z;
        			ge[6] = g[6]/z;
        			ge[7] = .000715487*y*y;
        			ge[8]=(-7.90744-.01792302*y+1.089687*z)*z;
        			ge[9] = -67.56741-2.46006*y;}
			    g[9] = Math.exp(g[9]);	// line 210
			    gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
			    if (kflag!=0) {
				    gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
			 	    gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}
			else {
	           if(z>3.4) error=true;
		       g[1] = -16.6904-.258318*y;
		       g[2] = (17.835+.154898*y)*z;
		       g[3] = (-.00971263+.0039774*z+.00009043*y)*y*y;
		       g[4] = (-5.94108-.0201335*y+.660432*z)*z*z;
		       g[5] = 85.469+11.7554*y;
		       g[6] = (-72.176-7.15723*y)*z;
		       g[7] = (-.041615+.0138147*z+.000545184*y)*y*y;
		       g[8] = (20.1758+1.0899*y-1.86438*z)*z*z;
		       g[9] = (288.3262+12.48536*y-88.16985*z-3.720309*y*z);
		       if ( kflag!=0) {
    			   gr[1]=-.258318;
	    		   gr[2] = .154898*z;
	        	   gr[3] = (-.01942526+.0079548*z+.00027129*y)*y;
		           gr[4] = -.0201335*z*z;
	    		   gr[5] = 11.7554;
	    		   gr[6] = -7.15723*z;
	     	       gr[7] = (-.08323+.0276294*z+.01635552*y)*y;
	    		   gr[8] = 1.0899*z*z;
	         	   gr[9] = 12.48536-3.720309*z;
   	    		   ge[2] = g[2]/z;
	    		   ge[3] = .0039774*y*y;
	    		   ge[4] = (-11.88216-.040267*y+1.981296*z)*z;
	    		   ge[6]= g[6]/z;
	    		   ge[7] = .0138147*y*y;
	    		   ge[8] = (40.3516+2.1798*y-5.59314*z)*z;
	    		   ge[9] = -88.16985-3.720309*y;}
			    g[9] = Math.exp(g[9]);	// line 210
			    gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
			    if (kflag!=0) {
				    gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
			 	    gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}}

		else if(z<=0.65) {gamm=1.396;if(kflag!=0) {gammr=0.0;gamme=0.0;}}
		else if(z<=1.5) {
    		g[1] = 1.39037-.00821107*y;
	    	g[2] = (-.0243213+.0124063*y)*z;
		    g[3] = (-.000725748-.00125272*z-.0000910474*y)*y*y;
    		g[4] = (.122447-.0189287*y-.148369*z)*z*z;
	    	g[5] = -3.75667-.154712*y;
		    g[6] = (7.78906+.203385*y)*z;
    		g[7] = (-.000716886+.00190409*z+.0000640147*y)*y*y;
	    	g[8] = (-5.34593-.0478112*y+1.27123*z)*z*z;
		    g[9] = Math.exp(86.0+5.058807*y-80.0*z-4.705882*y*z);
    		deno=1.0-g[9];
	    	gamm=g[1]+g[2]+g[3]+g[4]+(g[5]+g[6]+g[7]+g[8])/deno;
		    if (kflag!=0) {
    		    gr[1] = -.00821107;
    	    	gr[2] = .0124063*z;
	    	    gr[3] = (-.001451496-.00250544*z-.0002731422*y)*y;
    		    gr[4] = -.0189287*z*z;
    	    	gr[5] = -.154712;
	    	    gr[6] = .203385*z;
	        	gr[7] = (-.001433772+.00380818*z+.0001920441*y)*y;
        		gr[8] = -.0478112*z*z;
		        gr[9] = 5.058807-4.705882*z;
        		ge[2]=g[2]/z;
		        ge[3] = -.00125272*y*y;
        		ge[4] = (.244894-.0378754*y-.445107*z)*z;
		        ge[6]=g[6]/z;
        		ge[7] = .00190409*y*y;
		        ge[8] = (-10.69186-.0956224*y+3.81369*z)*z;
        		ge[9] = -80.0-4.705882*y;
		        gammr = gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/deno+(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/(deno*deno);
        		gamme = ge[2]+ge[3]+ge[4]+(ge[6]+ge[7]+ge[7]+ge[8])/deno+(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/(deno*deno);}}
        else if(z<=2.2) {
			g[1] = -17.0333-.508545*y;
			g[2] = (24.6299+.445617*y)*z;
			g[3] = (-.00895298+.00229618*z-.000289186*y)*y*y;
			g[4] = (-11.0204-.0989727*y+1.62903*z)*z*z;
			g[5] = 18.6797+.519662*y;
			g[6] = (-24.1338-.434837*y)*z;
			g[7] = (.00916089-.00152082*z+.000346482*y)*y*y;
			g[8] = (10.2035+.0970762*y-1.3946*z)*z*z;
			g[9] = (-142.762-1.647088*y+76.60312*z+.8259346*y*z);
			if (kflag!=0) {
				gr[1] = -.508545;
    			gr[2] = .445617*z;
	    		gr[3] = (-.01790596+.00459236*z-.000867558*y)*y;
		    	gr[4] = -.0989727*z*z;
			    gr[5] = .519662;
				gr[6] = -.434837*z;
    			gr[7] = (.01832178-.00304164*z+.001039446*y)*y;
	    		gr[8] = .0970762*z*z;
		    	gr[9] = -1.647088+.8259346*z;
				ge[2] = g[2]/z;
    			ge[3] = .00229618*y*y;
	    		ge[4] = (-22.0408-.1979454*y+4.88709*z)*z;
		    	ge[6] = g[6]/z;
				ge[7] = -.00152082*y*y;
    			ge[8] = (20.407+.1941524*y-4.1838*z)*z;
	    		ge[9] = 76.60312+.8259346*y;}
       	    g[9] = Math.exp(g[9]);	// line 210
	        gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
    	    if (kflag!=0) {
	   		    gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
	            gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}
        else if(z<=3.05) {
			g[1] = 2.24374+.103073*y;
    		g[2] = (-.532238-.0559852*y)*z;
	    	g[3] = (.00356484-.000101359*z+.000159127*y)*y*y;
		    g[4] = (-.0480156+.0106794*y+.0366035*z)*z*z;
			g[5] = -5.70378-.310056*y;
			g[6] = (5.01094+.180411*y)*z;
    		g[7] = (-.00949361+.00194839*z-.000224908*y)*y*y;
	    	g[8] = (-1.40331-.0279718*y+.120278*z)*z*z;
		    g[9] = (113.9755-4.985467*y-42.23833*z+2.009706*y*z);
			if (kflag!=0) {
				gr[1]=.103073;
        		gr[2] = -.0559852*z;
		        gr[3] = (.00712968-.00020218*z+.000477381*y)*y;
				gr[4] = .0106794*z*z;
        		gr[5] = -.310056;
		        gr[6] = .180411*z;
				gr[7] = (-.01898722+.00389678*z-.000674724*y)*y;
        		gr[8] = -.0279718*z*z;
		        gr[9] = -4.985467+2.009706*z;
				ge[2] = g[2]/z;
        		ge[3] = -.000101359*y*y;
		        ge[4] = (-.0960312+.0213588*y+.1098105*z)*z;
				ge[6] = g[6]/z;
        		ge[7] = .00194839*y*y;
		        ge[8] = (-2.80662-.0559436*y+.360834*z)*z;
				ge[9] = -42.23833+2.009706*y;}
    	    g[9] = Math.exp(g[9]);	// line 210
	        gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
		    if (kflag!=0) {
       		    gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
	            gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}
        else if(z<=3.4) {
			g[1] = -20.807+.40197*y;
        	g[2] = (22.591-0.2566*y)*z;
			g[3] = (-.00095833+.0023966*z+.00033671*y)*y*y;
    		g[4] = (-7.7174+.04606*y+0.878*z)*z*z;
    		g[5] = -217.37-4.6927*y;
        	g[6] = (181.01+2.6621*y)*z;
	        g[7] = (-.034759+.0064681*z-.00070391*y)*y*y;
			g[8] = (-50.019-.38381*y+4.5795*z)*z*z;
    		g[9] = (454.4373+12.50133*y-137.6001*z-3.641774*y*z);
	        if (kflag!=0) {
				gr[1] = 0.40197;
            	gr[2] = -0.2566*z;
				gr[3] = (-.00191666+.0047932*z+.00101013*y)*y;
            	gr[4] = .04606*z*z;
				gr[5] = -4.6927;
            	gr[6] = 2.6621*z;
				gr[7] = (-.069518+.0129362*z-.00211173*y)*y;
				gr[8] = -0.38381*z*z;
            	gr[9] = 12.50133-3.641774*z;
            	ge[2] = g[2]/z;
				ge[3] = .0023966*y*y;
            	ge[4] = (-15.4348+.09212*y+2.634*z)*z;
				ge[6] = g[6]/z;
            	ge[7] = .0064681*y*y;
				ge[8] = (-100.038-.76762*y+13.7385*z)*z;
            	ge[9] = -137.6001-3.641774*y;}
		    g[9] = Math.exp(g[9]);	// line 210
            gamm = g[1] + g[2] + g[3] + g[4] + (g[5] + g[6] + g[7] + g[8]) / (1.0+g[9]);
		    if (kflag!=0) {
    		    gammr=gr[1]+gr[2]+gr[3]+gr[4]+(gr[5]+gr[6]+gr[7]+gr[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*gr[9]*g[9]/((1.0+g[9])*(1.0+g[9]));
 	            gamme=ge[2] + ge[3]+ge[4]+(ge[6]+ge[7]+ge[8])/(1.0+g[9])-(g[5]+g[6]+g[7]+g[8])*ge[9]*g[9]/((1.0+g[9])*(1.0+g[9]));}}
        else {
            if(z>3.69) error=true;
			g[1] = -52.2951-.400011*y;
    		g[2] = (45.6439+.224484*y)*z;
	        g[3] = (-.00373775+.00243161*z+.000224755*y)*y*y;
			g[4] = (-12.9756-.0279517*y+1.22998*z)*z*z;
    		gamm=g[1]+g[2]+g[3]+g[4];
	        if (kflag!=0) {
				gr[1] = -.400011;
            	gr[2] = .224484*z;
				gr[3] = (-.0074755+.00486322*z+.000674265*y)*y;
            	gr[4] = -.0279517*z*z;
				ge[2] = g[2]/z;
            	ge[3] = .00243161*y*y;
				ge[4] = (-25.9512-.0559034*y+3.68994*z)*z;
            	gammr = gr[1] +gr[2]+gr[3]+gr[4];
				gamme = ge[2]+ge[3]+ge[4];}}

		p = (gamm-1)*e*r;	// line 230
		if (kflag!=0) {
			gammr/=2.302585;
			gamme/=2.302585;
			a = e*((gamm-1.0)*(gamm+gamme)+gammr);}

        if(iflag==-1 || jflag==1) return(0);
        if(jflag==0) {
            plow=p;alow=a;
			p = plow+(phigh-plow)/(yhigh-ylow)*(ym-ylow);
			a = alow+(ahigh-alow)/(yhigh-ylow)*(ym-ylow);
            r=rsave;
            return(0);}
        phigh=p;ahigh=a;y=-4.5-0.025;
        if(iflag==1) y=-.5-.05;
        ylow=y;r=Math.pow(10.,y)*r0;
        jflag=0;}
    }

int t_pr() {	// temperature as a function of pressure and density (old tgas3.f)

	double r0 = 1.292,p0 = 101330,t0 = 273.15,gascon = 287.06;
	double y,x,tlow,thigh=0.,ylow=0.,yhigh=0.,ym=0.,rsave=0.,rho,z1,tnon;
	int iflag,jflag=0;

    rho=r;
	y=(Math.log(rho/r0))/(Math.log(10));
	x=(Math.log(p/p0))/(Math.log(10));


	if (Math.abs(y+4.5)<.045){iflag=0;rsave= rho;ym = y;y= -4.5+.045;yhigh = y;rho = Math.pow(10.0,y)*r0;jflag = -1;}
	else if (Math.abs(y+0.5)<.005){iflag = 1;rsave = rho;ym = y;y = -0.5 + .005;yhigh =y;rho = Math.pow(10.0,y)*r0;jflag=-1;}
	else iflag = -1;

	for(;;) {
		z1 = x-y;
		if (y>-0.5){			// goto 190
        	if (z1<= .25) {
        		t = p/(rho * gascon);}
        	else if(z1 <= 1) {
        		tnon = (-.00154141+.000658337*y) + ((.982201-.00385028*y)*z1) + ((.000123111-.00040821*z1+.0000213592*y)*y*y) + ((.0377441+.00456963*y-.0235172*z1)*z1*z1);
        		t = Math.pow(10,tnon)*t0;}
        	else if(z1<=1.45) {
        		tnon = (.806492+.0991293*y) + ((-1.70742-.228264*y)*z1) + ((.005035-.00613927*z1+.000169824*y)*y*y) + ((3.02351+.131574*y-1.12755*z1)*z1*z1) + ((-.11793-.212207*y) + ((1.36524+.405886*y)*z1) + ((-.018823+.0165486*z1-.0005114*y)*y*y) + ((-2.10926-.189881*y+.879806*z1)*z1*z1))/(1.0+(Math.exp(195.9604-42.69391*y-173.4931*z1+37.62898*y*z1)));
        		t = Math.pow(10, tnon)*t0;}
        	else {
        	    if(z1>2.3) error=true;
        		tnon = (-1.66249-.0891113*y) + ((4.11648+.0878093*y)*z1) + ((-.00309742+.00199879*z1+.0000685472*y)*y*y) + ((-1.84445-.00750324*y+.305784*z1)*z1*z1) + ((11.1555+1.321*y) + ((-17.1236-1.2919*y)*z1) + ((.0628124-.0307949*z1+.00157743*y)*y*y) + ((8.63804+.307809*y-1.42634*z1)*z1*z1))/(1.0+(Math.exp(133.0611+8.979635*y-72.65298*z1-2.449009*y*z1)));
        		t = Math.pow(10, tnon)*t0;}}
		else if (y > -4.5) {  			// goto 140
        	if (z1 <= 0.25) {
        		t = p / (rho*gascon);}
        	else if( z1 <= 0.95) {
        		tnon = (.020391+.0076731*y) + ((.848581-.0293086*y)*z1) + ((.000840269-.00147701*z1+.0000313687*y)*y*y) + ((0.267251+.0237262*y-.141973*z1)*z1*z1);
        		t = (Math.pow(10.0,tnon)*t0);}
        	else if (z1 <= 1.45) {
        		tnon = (-5.12404-.2847*y) + ((15.4532+.452475*y)*z1) + ((-.0122881+.00856845*z1-.000325256*y)*y*y) + ((-13.5181-.168725*y+4.18451*z1)*z1*z1) + ((7.52564+.835238*y) + ((-19.5558-1.23393*y)*z1) + ((.033451-.0234269*z1+.000481788*y)*y*y) + ((17.1779+.454628*y-5.09936*z1)*z1*z1))/(1.0+(Math.exp(61.48442-18.28123*y-54.68755*z1+15.625*y*z1)));
        		t = Math.pow(10, tnon)*t0;}
        	else if (z1 <= 2.05) {
        		tnon = (-12.3779-1.14728*y) + ((24.1382+1.38957*y)*z1) + ((-.0363693+.0224265*z1-.000323888*y)*y*y) + ((-14.2844-.406553*y+2.8752*z1)*z1*z1) + ((4.40782+1.33046*y) + ((-11.545-1.59892*y)*z1) + ((.053058-.0310376*z1+.00047765*y)*y*y) + ((8.57309+.471274*y-1.96233*z1)*z1*z1))/(1.0+(Math.exp(140.75-6.499992*y-77.5*z1+5.0*y*z1)));
        		t = Math.pow(10, tnon)*t0;}
        	else {
        	    if (z1>2.5) error=true;
        		tnon = (-12.7244-1.66684*y) + ((17.2708+1.45307*y)*z1) + ((-.0364515+.0190463*z1+.000480787*y)*y*y) + ((-6.97208-.304323*y+.967524*z1)*z1*z1) + ((7.7133+.50834*y) + ((-9.8211-.449138*y)*z1) + ((-.000941787-.00240293*z1-.00082845*y)*y*y) + ((4.1653+.0963923*y-.588807*z1)*z1*z1))/(1.0+(Math.exp(-1092.654-305.312*y+465.6243*z1+131.2498*y*z1)));
        		t = Math.pow(10, tnon)*t0;}}
		else if (z1 > 0.25) { 			 // goto 100
        	if (z1 <= 0.95) {
        		tnon = (.123718+.0108623*y) + ((.224239-.0824608*y)*z1) + ((-.00117615-.00187566*z1-.000119155*y)*y*y) +((1.18397+.0648520*y-.55263*z1)*z1*z1);
        		t = Math.pow(10.0,tnon) * t0;}
        	else if (z1 <= 1.4) {
        		tnon = (-8.12952-.828637*y) + ((22.6904+1.41132*y)*z1) + ((-.0298633+.0270066*z1-.000228103*y)*y*y) + ((-19.1806-.578875*y+5.6258*z1)*z1*z1) + ((-3.99845+.226369*y) + ((2.52876-.728448*y)*z1) + ((.0109769-.0183819*z1-.00015138*y)*y*y) + ((2.99238+.39144*y-2.04463*z1)*z1*z1))/(1.0+(Math.exp(-38.87015-29.08228*y+40.70557*z1+26.82347*y*z1)));
        		t = Math.pow(10.0,tnon) * t0;}
        	else if (z1 <= 1.95) {
        		tnon = (-19.8573-1.67225*y) + ((37.6159+2.10964*y)*z1) + ((-.0340174+.0231717*z1-.0000980257*y)*y*y) + ((-22.2215-.644596*y+4.40486*z1)*z1*z1) + ((-5.36809+.241201*y) + ((-1.25881-.862744*y)*z1) + ((-.00379774-.00781335*z1-.000380005*y)*y*y) + ((5.58609+.378963*y-1.81566*z1)*z1*z1))/(1.0+(Math.exp(20.8-25.6*y+1.0*z1+18.0*y*z1)));
        		t = Math.pow(10.0,tnon) * t0;}
        	else {
        	    if (z1>2.6) error=true;
        		tnon = (-23.3271-1.89958*y) + ((32.144+1.68622*y)*z1) + ((-.0442123+.0282629*z1+.000663272*y)*y*y) + ((-13.8645-.340976*y+2.04466*z1)*z1*z1) + ((8.35474+1.71347*y) + ((-16.0715-1.63139*y)*z1) + ((.0414641-.0230068*z1+.0000153246*y)*y*y) + ((8.70275+.360966*y-1.46166*z1)*z1*z1))/(1.0+(Math.exp(111.5884-6.452606*y-53.37863*z1+2.026986*y*z1)));
        		t = Math.pow(10.0,tnon) * t0;}}
		else t = p / (rho * gascon);

        if(iflag==-1 || jflag==1) return(0);
        if(jflag==0) {
			tlow = t;
			t = tlow + (thigh - tlow) / (yhigh - ylow) * (ym - ylow);
			rho = rsave;
			return(0);}
		thigh = t;y = -4.5 - .045;
		if(iflag==1) y = -0.5 - 0.005;
		ylow = y;
		rho = Math.pow(10,y) * rho;
	    jflag = 0;}
	}

}