s_re.m

function[s]=s_re(r,e)
% entropy as a function of energy and density (old tgas2.f)

e0 = 78408.4;
r0 = 1.292;
gascon = 287.06;
eratio = e/e0;
rratio = r/r0;
y = log(rratio)/log(10);
z = log(eratio)/log(10);
error=false;

if(abs(y+4.5)<.045) 
    iflag=0;
    jflag=-1;
    rsave=r;
    ym=y;
    y=-4.5+.045;
    yhigh=y;
    r=(10^y)*r0;
elseif(abs(y+0.5)<.005)
    iflag=1;
    jflag=-1;
    rsave=r;
    ym=y;
    y=-0.5+0.005;
    yhigh=y;
    r=(10^y)*r0;
else
    iflag=-1;
end

while(1)
	  	if (z<=0.65) 
			deltz = z-0.4;
			delts=(2.5*deltz-y)*gascon*2.302585;
			s = 6779.2004+delts;        
		elseif(y>-4.5) 
			if(y>-0.5) 
			    if(z>3.0) 
                error=true;
                end
			    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;
                end
		    	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));
            end
        else
            if(z>3.69) 
                error=true;
            end
		    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));
        end

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