r_pt.m

function[r]=r_pt(p,t)         
% Determines density iteratively as functions of p and t (old airtbn routine)     
    
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:1:10
        %Bracket solution;
            for j=1:1:2
                r=rr(j)*r0;
                t=t_pr(p,r);
                gg(j)=(t-temp)/temp;
            end

            if(gg(1)*gg(2)>0.0) 
                rr(1)=rr(1)*2.;
                rr(2)=rr(2)*0.5;
            else
                break;
            end
    end 
            
          if(ic>=10) 
              error=true;
              return;
          end

          for j=1:1:50
          %Regular Falsi
            rrn=(rr(2)*gg(1)-rr(1)*gg(2))/(gg(1)-gg(2));
            r=rrn*r0;
            t=t_pr(p,r);
            ggn=(t-temp)/temp;
            aggn=abs(ggn);

            if(aggn<0.1) 
                break
            end 

            if(ggn*gg(1)>0.0) 
                rr(1)=rrn;
                gg(1)=ggn;
            else
                rr(2)=rrn;
                gg(2)=ggn;
            end
          end                 
          
          for j=3:1:50
          % 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=t_pr(p,r);
            gg(j)=(t-temp)/temp;

            if(abs(gg(j))<5.e-4) 
                break;
            end
          end
          
          t=temp;
          error=false;
          return;
       
end