//************
//  Description: Create a data ramp in a Midas file
//
//  Author:	JGS   4//90
//
//  $Id: ramp.cc,v 1.1 2009/08/11 23:42:59 ntn Exp $
//  
//************

#include <limits>
#include <stdexcept>
#include <vector>
#include <primitive.h>

void mainroutine () 
{
  CPHEADER HCB;

  // test for allowable formats (any scalar numeric)
  string filename = m_apick(1);
  string format = m_test_format(m_apick(2), "S#,-SP");

  // fill out header and open ouput file
  m_init(HCB, filename, "1000", format, 0);
  HCB.xstart = m_dpick(6);
  HCB.xdelta = m_dpick(7);
  HCB.xunits = m_bpick(8);
  HCB.size = m_dpick(3);
  m_allocate(HCB);

  // get start/stride 
  real_8 dstart = m_dpick(4);
  real_8 ddelta = m_dpick(5);

  // find buffer size
  size_t LPER = m_xfer_length(HCB.bpe, -8192);
  vector<int_1> buf(m_bytes_reqd(HCB.bpe, LPER));
  HCB.buf_type = m_promote_type(HCB.buf_type, 'B');

  // signal start of processing block
  m_sync();

  // loop on buffers depending on data type
  m_todo(HCB.size);
  int_4 NN;
  while (m_do(LPER, NN)) {
    
    switch (HCB.buf_type) {
    case 'D': {
	real_8* dbuf = reinterpret_cast<real_8*>(&buf[0]);
	for (int ii = 0; ii < NN; ++ii) {
	  dbuf[ii] = dstart;
	  dstart += ddelta;
	}
      }
      break;
    case 'F': {
	real_4* fbuf = reinterpret_cast<real_4*>(&buf[0]);
	for (int ii = 0; ii < NN; ++ii) {
	  fbuf[ii] = dstart;
	  dstart += ddelta;
	}
      }
      break;

      // ///// Integer types

      // The integer output of RAMP is presently
      // implementation-dependent.  The standard only defines
      // floating-to-integral and integral-to-integral conversions
      // when the target integer can represent the source value. For
      // many years, RAMP has depended on a common implementation
      // detail of integral-to-integral conversions: wraparound
      // arithmetic, so that assigning an int_4 like 32770 to an int_2
      // produces -32766.

    case 'X': {
	int_8* xbuf = reinterpret_cast<int_8*>(&buf[0]);
	for (int ii = 0; ii < NN; ++ii) {
	  xbuf[ii] = int_8(dstart > 0  ?  dstart + 0.5  :  dstart - 0.5);
	  dstart += ddelta;
	}
      }
      break;
    case 'L': {
	int_4* lbuf = reinterpret_cast<int_4*>(&buf[0]);
	for (int ii = 0; ii < NN; ++ii) {
	  lbuf[ii] = int_8(dstart > 0  ?  dstart + 0.5  :  dstart - 0.5);
	  dstart += ddelta;
	}
      }
      break;
    case 'I': {
	int_2* ibuf = reinterpret_cast<int_2*>(&buf[0]);
	for (int ii = 0; ii < NN; ++ii) {
	  ibuf[ii] = int_4(dstart > 0  ?  dstart + 0.5  :  dstart - 0.5);
	  dstart += ddelta;
	}
      }
      break;
    case 'B':
      // Although technically an int_2 is all that's required here, we
      // use the int_4 cast, because it's very possible to have an
      // internal buffer that spans more than 65536. int_4 goes up to
      // 2G (signed), which is better.
      for (int ii = 0; ii < NN; ++ii) {
        buf[ii] = int_4(dstart > 0  ?  dstart + 0.5  :  dstart - 0.5);
        dstart += ddelta;
      }
      break;
    default:
      throw logic_error(string("No handling for buffer type ") +
			HCB.buf_type);
    } // switch(HCB.buf_type)
    
    m_filad(HCB, &buf[0], NN);

    // For integer types, don't waste mantissa bits on range that
    // can't be represented in the output type. Bring dstart back
    // into range after each call to m_filad.

    // Note that in the code below, we test for dstart being outside
    // of the signed range, but need to correct it by the full range,
    // which is the max value of the unsigned equivalent, plus one
    // (last - start + 1 == size).

    switch (HCB.buf_type) {
    case 'L':
      while (dstart > numeric_limits<int_4>::max())
	dstart -= numeric_limits<int_u4>::max() + 1.0;
      while (dstart < numeric_limits<int_4>::min())
	dstart += numeric_limits<int_u4>::max() + 1.0;
      break;
    case 'I':
      while (dstart > numeric_limits<int_2>::max())
	dstart -= numeric_limits<int_u2>::max() + 1.0;
      while (dstart < numeric_limits<int_2>::min())
	dstart += numeric_limits<int_u2>::max() + 1.0;
      break;
    case 'B':
      while (dstart > numeric_limits<int_1>::max())
	dstart -= numeric_limits<int_u1>::max() + 1.0;
      while (dstart < numeric_limits<int_1>::min())
	dstart += numeric_limits<int_u1>::max() + 1.0;
      break;
    } // switch(HCB.buf_type)

  }

  // close output file
  m_close(HCB);
}

