function [Results,keep,etaResults,gResults,etaLResults]=ideamodelsolution(Omega,Beta,rho,PhiH,PhiL,Gamma,SigmaH,SigmaL,theta,lambda,ZBar,Lbar,MBar,DltaBar,Tau,PBar,qbar);

% function [Results,keep]=ideamodelsolution(Omega,Beta,rho,PhiH,PhiL,Gamma,SigmaH,SigmaL,theta,lambda,ZBar,Lbar,MBar,DltaBar,Tau,PBar,qbar);
%
%  Basic function called by IdeaModel*.m (and others) to solve
%  the model for various combinations of parameter values
%
%  6/17/15 -- Two state Luttmer/GLLM model with muH, muL, sigmaH, sigmaL, pbar, qbar
%             (and no reflecting barrier so Ex and Elogx are easy to compute)

NumPeriods=length(PhiH);
gResults=zeros(NumPeriods,1)*NaN;
etaResults=gResults; sResults=gResults; dResults=gResults; muHResults=gResults; etaLResults=gResults;
eResults=zeros(NumPeriods,2)*NaN;

                       
sTries=[.1 .05 .15 .02 .2 .3];  % A range of initial guesses

for t=1:NumPeriods;
    f=@(ss) e_solveideamodel(ss,Omega,Beta(t),rho,PhiH(t),PhiL(t),Gamma(t),SigmaH(t),SigmaL(t),theta,lambda,ZBar(t),Lbar(t),MBar(t),DltaBar(t),Tau(t),PBar(t),qbar);
    fresults=@(ss) solveideamodel(ss,Omega,Beta(t),rho,PhiH(t),PhiL(t),Gamma(t),SigmaH(t),SigmaL(t),theta,lambda,ZBar(t),Lbar(t),MBar(t),DltaBar(t),Tau(t),PBar(t),qbar);
    %tresults=@(ss) testsolvemodel(ss,Omega,Beta(t),rho,PhiH(t),PhiL(t),Gamma(t),SigmaH(t),SigmaL(t),theta,lambda,ZBar(t),Lbar(t),MBar(t),DltaBar(t),Tau(t),PBar(t),qbar);

    % Ensure conditions are all satisfied:
    %
    %   1. mu>0 and 0<s<1
    %   2. f(s-eps)>0 and f(s+eps)<0  -- this picks out the "right" solution(?)
    keeptrying=1; i=1;
    while keeptrying;
        %fprintf('Try #%1f: sguess=%4.2f\n',[i sTries(i)]);
        if isreal(f(sTries(i)));
            s=fzero(f,sTries(i));
            [sblah estarH estarL g eta dlta muH muL]=fresults(s);
            if s<1 & s>0 & muH>0;
                eleft=f(s-.0001);
                eright=f(s+.0001);
                if eleft>0 & eright<0; keeptrying=0; end; 
            end; 
        end; 
        i=i+1;
        if i>length(sTries); disp '****** No solution found... Setting s=1.05 ******';
            s=1.05; % feed something to flag the error break 
                    %testfzero(f,(0.07:.01:0.09));
                    %keyboard; 
            break
        end; 
    end;
    [sblah,estarH,estarL,g,eta,dlta,muH,muL]=solveideamodel(s,Omega,Beta(t),rho,PhiH(t),PhiL(t),Gamma(t),SigmaH(t),SigmaL(t),theta,lambda,ZBar(t),Lbar(t),MBar(t),DltaBar(t),Tau(t),PBar(t),qbar);

    cH=Beta(t)*log(Omega-estarH)+(muH-.5*SigmaH(t)^2)/(rho+dlta+DltaBar(t));
    cL=Beta(t)*log(Omega-estarL)+(muL-.5*SigmaL(t)^2)/(rho+dlta+DltaBar(t));
    if cH<cL;
        disp ' ';
        fprintf('Problem in ideamodelsolution.m:  cH = %8.4f    cL = %8.4f\n',[cH cL]);
    end;

    % if mu<0; p=-99; end; % Throw out if estar<0...
    eResults(t,:)=[estarH estarL];
    gResults(t)=g*100;
    etaResults(t)=eta;
    sResults(t)=s;
    dResults(t)=dlta; 
    muHResults(t)=muH;
    fprintf('muL = %8.4f, estarL = %8.4f      cH = %8.4f, cL = %8.4f\n',[muL estarL cH cL]);
    
    % Check Low condition
    mutildeL=muL-.5*SigmaL(t)^2;
    mm=mutildeL/SigmaL(t)^2;
    etaL = 1/(-mm+sqrt(mm^2+2*DltaBar(t)/SigmaL(t)^2));
    etaLResults(t)=etaL;
    if etaL>eta; fprintf('Wait! etaL=%6.4f > eta=%6.4f for these parameter values!!!!\n',[etaL eta]); end;
    
    % %%% testfzero(f,(0.06:.001:.105))    
    %paramrange=(0.00:.0001:.21);
    %test1
    %wait;
    %keyboard
    
end; % t


% Now, show all collections that feature 0<eta<1 -- see "linear indexing" in matlab
disp ' ';
mutildeH=muHResults-1/2*SigmaH.^2;
indx=find(etaResults>0);
fmt='%6.0f %7.3f';
tle='t eta gx100 mutldeH s beta gamma zbar phiH phiL sigmaH dlta dltabar tau PBar eH eL MBar Lbar';
data=[1979+(1:NumPeriods)' etaResults gResults mutildeH sResults Beta Gamma ZBar PhiH PhiL SigmaH dResults DltaBar Tau PBar eResults MBar Lbar];
cshow(' ',data,fmt,tle);
s=data(:,4);
keep=s~=1.05;

Results=data;
