/**
    @file MergeRefrac.cpp
    @author Brian Magill
    @datecreated 08/13/2008

    $Date:$
    $Revision:$

    @copyright
     (©) Copyright 2006 by GATS Inc.
     11864 Canon Blvd., Suite 101, Newport News, VA 23606

     All Rights Reserved. No part of this software or publication may be
     reproduced, stored in a retrieval system, or transmitted, in any form
     or by any means, electronic, mechanical, photocopying, recording, or
     otherwise without the prior written permission of GATS Inc.

    @brief: Merges the two refraction profiles  
    
*/
//----------------------------------------------------------------------
//

#include "MergeRefrac.h"
#include <string>
#include <stdexcept>
#include <cassert>
#include <iostream>

using namespace std;

long MergeRefrac::FindIndex(const valarray<double> & array) const
{
long index = 0;
string strErr;

if (mergeCriteria < array.min() || array.max() <  mergeCriteria)
    {
    strErr = "MergeRefrac::findIndex error: mergeCriteria out of range";
    throw runtime_error(strErr);
    }    
for (long i = 0; i < array.size(); i++)
    {
    if (mergeCriteria < array[i])
        {
        index = i;  
        break;  
        }
    }
    
return index;

}


valarray<double> MergeRefrac::Combine(const valarray<double> &ProfileA, 
                                            const valarray<double> &ProfileB) const
{
assert(ProfileA.size() == ProfileB.size());
long profileSize = ProfileA.size();

valarray<double> work(profileSize);
valarray<double> weightB(0.0, profileSize);    
long i;

long index = FindIndex(ProfileB);

if (index - halfWidthPoints < 1 || index + halfWidthPoints > profileSize)
    throw runtime_error("MergeRefrac::Combine error: merge interval out of range");

double m = 2.0*(halfWidthPoints + 1);
//long offset = index - halfWidthPoints;

for(i = 0; i < 2*halfWidthPoints + 1; i++)
    { weightB[i + index - halfWidthPoints - 1] = (i + 1)/m; }
    
for(i = index + halfWidthPoints; i < weightB.size(); i++)
    { weightB[i] = 1.0; }


valarray<double> weightA = valarray<double> (1.0, profileSize) - weightB;    

// cout << endl;
// for(i = 0; i < weightB.size(); i++)
//     { cout << ProfileA[i] << "\t" << ProfileB[i] << "\t" << weightA[i] << "\t" << weightB[i] << endl; }

work = weightA*ProfileA + weightB*ProfileB;
   
return work;
}