function [obs] = mcmc_fwd(np,no,params,obs_loc,model)
%...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW
%----------------------------------------------------------------------------------85
%                            ~ ~ ~  MCMC_FWD.M ~ ~ ~                                %
%-----------------------------------------------------------------------------------%
%...DESCRIPTION:                                                                    %
%...  Matlab function which computes a vector of forward observations using some    %
%...  user-specified (and coded) forward model and values of model control          %
%...  parameters.                                                                   %
%...                                                                                %
%...INPUTS:                                                                         %
%...  np: Number of parameters                                                      %
%...  no: Number of observations                                                    %
%...  params: Array of parameter values                                             %
%...  obs_loc: Array specifying 'location' of observation (usually a spatial gridpt)%
%...  model: Forward model whose parameters will be inverted by MCMC                %
%...      'sine'    ......   Sinusoid                                               %
%...      '1dadv'   ......   1D advection - calls advection, RK4adv                 %
%...                                                                                %
%...OUTPUT:                                                                         %
%...  obs: Array of forward observation values.                                     %
%...                                                                                %
%...UPDATES:                                                                        %
%...  April 12 2011     ... work started                                            %
%-----------------------------------------------------------------------------------%
%...Created by Marcus van Lier-Walqui for NSF-supported work (PI: Sharan Majumdar & %
%...Tomislava Vukicevic).                                                           %
%...RSMAS, University of  Miami, Miami FL.                                          %
%...                                                                                %
%...COME ON YOU GUNNERS!!!                                                          %
%...                                                                                %
%-----------------------------------------------------------------------------------%
%...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW...MvLW





%------------------------------------ 1D Advection ---------------------------------%
if strcmp(model,'1dadv')
  %----------------- Initialize Model Parameters -----------------%
  %...Note that some of these definitions will be redundant with those defined within
  %...The main code
  obs = zeros(1,no);  %...'allocate' observational array
  
  %--------------------- Advection Model Settings ----------------%
  
  ntime   = 500;            %...Number of time steps
  dt      = 10.;            %...Time step
  nspace  = 100;            %...Number of spatial grid intervals
  L       = 36000.;         %...Domain size (in some distance units)
  
  
  ds=L/nspace;              %...Spatial grid step
  
  %...Spatial grid location values
  sgrid = zeros(1,nspace+1);
  for ns=2:nspace
    sgrid(ns)=(ns-1)*ds;
  end
  sgrid(nspace+1)=L;
  
  %...Set initial condition parameters for two-wave state superimposed on a basic 
  %...state. 
  %...          u_init = a0 + a1*sin(wn1*s) + a2*cos(wn2*s)
  %...
  %...Wavelenths
  [wl]=[params(1),params(2)]*L;
  %...Wavenumbers
  [wn]=[2*pi/wl(1,1),2*pi/wl(1,2)];
  %...Wave amplitudes
  [A]=[params(4),params(5)];
  %...Basic state amplitude
  A0=params(3);

  %...Compute initial condition
  u0 = zeros(1,nspace+1);
  for ns=1:nspace
    u0(ns)=A0+A(1,1)*sin(wn(1,1)*sgrid(ns))+A(1,2)*cos(wn(1,2)*sgrid(ns)); 
  end
  u0(nspace+1)=u0(1);
  %.................................................................................%
  %...                END OF `USER INPUT'
  %.................................................................................%
  %...                BEGIN MODEL SIMULATION
  %.................................................................................%
  %...
  [up]=advection(u0,nspace,dt,ds,ntime);
  %...
  %...Add periodic point (?)
  for i=1:nspace+1
    up(1,i)=u0(1,i);
  end

  %...Fill forward observation vector
  for i=1:no
    obs(i) = up(obs_loc(i,1),obs_loc(i,2));
  end

  %obs(1) = up(obs_loc(1,1),obs_loc(1,2));
  %obs(2) = up(obs_loc(2,1),obs_loc(2,2));
  %obs(3) = up(obs_loc(3,1),obs_loc(3,2));

%------------------------------------ Gaussian -------------------------------------%
elseif strcmp(model,'gaussian')
  %...All this does is return the same value of the parameter as observation. 
  obs = zeros(1,no);
  for i=1:np
    for j=1:no
      obs(j) = obs(j)+obs_loc(j,i)*(params(i));
    end
  end


%--------------------------- Damped Harmonic Oscillator ----------------------------%
%...Must allow for either initial amplitude or coefficient inversion.
elseif strcmp(model,'dampharm')
  %...Harmonic oscillator settings
  tspan = 0:1:50;
  init  = [params(1),0];

  pprr  = @(t, y) mcmc_dampharm(t,y,params(2),params(3));
  [T,Y] = ode45(pprr,tspan,init);

  for i=1:no
    obs(i) = Y(obs_loc(i));
  end

%--------------------------- Lorenz 3-Component Model (1963) -----------------------%
elseif strcmp(model,'lorenz')
  %...Lorenz model settings
  tspan = 0:.02:4.98;
  init  = [params(1) params(2) params(3)];
  p = [params(4) params(5) params(6)];

  [t,x] = ode45(@mcmc_lorenz,tspan,init,[],p);
  for i=1:no
    obs(i)  = x(obs_loc(i,1),obs_loc(i,2));
  end
  

else
  'DOES NOT COMPUTE!'
  return
end