% ionomodel.m: user-supplied ionospheric model
% GUISDAP v.1.60 96-05-27 Copyright Asko Huuskonen and Markku Lehtinen
% 
% 'ionomodel' is a user-supplied function which for a given set of heights
% outputs the plasma parameters (Ne, Ti, Te/Ti, coll, [O+]/Ne, velocity)
% with their a priori uncertainties. This simple example is static, a more complete
% model might use the date and time values to modify the model.
%
% function [apriori, apriorierror]=ionomodel(heights);

function [apriori, apriorierror]=ionomodel(heights);

global path_GUP ionomodel_first ionomodel_Ne_mult

%maxCollAltit=107; % Below this collision frequency is given freedom to change
%minRatioAltit=107; % Above this temperture ratio is given freedom to change
maxCollAltit=100 % Below this collision frequency is given freedom to change
minRatioAltit=100; % Above this temperture ratio is given freedom to change
% Note: This above need not be equal, but setting minRatioAltit less than
% maxCollAltit will produce unconvergent fits, because both parameters will change.
% It is completely OK to have a range where neither is fitted.

%--- Appended by S. Saito ---
global MCA MRA
MCA=maxCollAltit;
MRA=minRatioAltit;
%----------------------------

len=length(heights);
heights=col(heights);

if nargout==2, error=1; else, error=0; end

file=canon([path_GUP,'matfiles/summer1996'],0);
load(file)
if isempty(ionomodel_first)
  fprintf('*\n*\n* This ionomodel uses the GIVEME model results for Summer 1996\n*\n*\n')
  ionomodel_first=1;
end
heights2=heights;
heights=max(heights,altitude(1));
heights=min(heights,altitude(length(altitude)));

if ~isempty(ionomodel_Ne_mult)
  Nemult=ionomodel_Ne_mult;
else
  Nemult=1;
end

par=1; % Electron density
   apriori(:,par)=Nemult*exp(inter3(heights2,altitude,log(n_electr)))';
   if error, apriorierror(:,par)=1000*apriori(:,par); end

c_O=1-(n_NO_ion+n_02_ion+n_N2_ion)./n_electr;
par=2; % Ion temperature
   apriori(:,par)=inter3(heights,altitude,c_O.*t_O__ion+(1-c_O).*t_M2_ion)';
   if error, apriorierror(:,par)=1000*apriori(:,par); end
  
par=3; % Temperature ratio
   te=inter3(heights,altitude,t_electr)';
   apriori(:,par)=te./apriori(:,2);
   if error,
    apriorierror(:,par)=1000*apriori(:,par);
    ind=find(heights<=minRatioAltit);
    if length(ind)>0,apriorierror(ind,par)=apriori(ind,par)/1000;end
  end
  
par=4; % Ion-neutral collision frequency
   apriori(:,par)=max(3578*(exp(-(heights2-100)/5.8)),10);
   if error,
     apriorierror(:,par)=apriori(:,par)/100;
     ind=find(heights<=maxCollAltit);
     if length(ind)>0,apriorierror(ind,par)=apriori(ind,par);end
   end
 
par=5; % Ion velocity
  apriori(:,par)=zeros(len,1);
  if error, apriorierror(:,par)=100000*ones(len,1); end

par=6; % Ion composition
   apriori(:,par)=inter3(heights',altitude,c_O)';
   if error, apriorierror(:,par)=ones(len,1); end

