#include "m_pd.h"
#include "math.h"

#define UNITBIT32 1572864.  /* 3*2^19; bit 32 has place value 1 */
/*    should have been by CPU type and not by operating system! */
#ifdef IRIX
    /* big-endian.  Most significant byte is at low address in memory */
#define HIOFFSET 0    /* word offset to find MSB */
#define LOWOFFSET 1    /* word offset to find LSB */
#define int32 long  /* a data type that has 32 bits */
#else
#ifdef NT
    /* little-endian; most significant byte is at highest address */
#define HIOFFSET 1
#define LOWOFFSET 0
#define int32 long
#else
#ifdef __linux__

#include <endian.h>

#if !defined(__BYTE_ORDER) || !defined(__LITTLE_ENDIAN)                         
#error No byte order defined                                                    
#endif                                                                          
                                                                                
#if __BYTE_ORDER == __LITTLE_ENDIAN                                             
#define HIOFFSET 1                                                              
#define LOWOFFSET 0                                                             
#else                                                                           
#define HIOFFSET 0    /* word offset to find MSB */                             
#define LOWOFFSET 1    /* word offset to find LSB */                            
#endif /* __BYTE_ORDER */                                                       

#include <sys/types.h>
#define int32 int32_t

#else
#ifdef MACOSX
#define HIOFFSET 0    /* word offset to find MSB */
#define LOWOFFSET 1    /* word offset to find LSB */
#define int32 int  /* a data type that has 32 bits */

#endif /* MACOSX */
#endif /* __linux__ */
#endif /* NT */
#endif /* SGI */



union tabfudge
{
    double tf_d;
    int32 tf_i[2];
};



static t_class *buzz_tilde_class;


typedef struct _buzz_tilde
{
    t_object x_obj;
    double x_phase;
    float x_conv;
    t_float x_a;
    t_float x_H;
    float _S1;
    float _S2;
    float x_f;	    /* scalar frequency */
} t_buzz_tilde;

static void *buzz_tilde_new(t_floatarg f)
{
    t_buzz_tilde *x = (t_buzz_tilde *)pd_new(buzz_tilde_class);
    x->x_f = f;
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_float, gensym("ft1"));
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_float, gensym("f_a"));
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_float, gensym("f_H"));
    x->x_phase = 0;
    x->x_conv = 0;
    x->x_a = 0;
    x->x_H = 0;
    x->_S1 = 0;
    x->_S2 = 0;

    outlet_new(&x->x_obj, gensym("signal"));
    return (x);
}

static t_int *buzz_tilde_perform(t_int *w)
{
    t_buzz_tilde *x = (t_buzz_tilde *)(w[1]);
    t_float *in = (t_float *)(w[2]);
    t_float *out = (t_float *)(w[3]);
    int n = (int)(w[4]);
    double dphase = x->x_phase + UNITBIT32;
    union tabfudge tf;
    int normhipart;
    float conv = x->x_conv;
    
    double my_phi;
    float a = x->x_a;
    float S1 = x->_S1;
    float S2 = x->_S2;
    float N = x->x_H + 2;
    float pi2 = 2.f * M_PI;
    float a2 = 2 * a ;
    float ap2 = a*a;
	   
    tf.tf_d = UNITBIT32;
    normhipart = tf.tf_i[HIOFFSET];
    tf.tf_d = dphase;

    while (n--)
    {
    	tf.tf_i[HIOFFSET] = normhipart;
    	// increment phase
	dphase += *in++ * conv;
	my_phi = pi2 * dphase;
	
	// do work with my_phi:
	*out++ = ( 
			( 1 - a*cos(my_phi) ) * 
			( S1*cos(my_phi) - S2*cos(N*my_phi) ) -
			a * sin(my_phi) * 
			(  S2*sin(N*my_phi) - S1 * sin(my_phi) 	)
		) 
		/ 
		( 1 - a2 * cos(my_phi) + ap2 ) ;
		
	// DON'T: send old phase to signal outlet
	//*out++ = tf.tf_d - UNITBIT32;
	tf.tf_d = dphase;
    }
    tf.tf_i[HIOFFSET] = normhipart;
    x->x_phase = tf.tf_d - UNITBIT32;
    return (w+5);
}

static void buzz_tilde_dsp(t_buzz_tilde *x, t_signal **sp)
{
    x->x_conv = 1./sp[0]->s_sr;
    dsp_add(buzz_tilde_perform, 4, x, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
}

static void varsetup (t_buzz_tilde *x)
{
    float a = x->x_a;
    float H = x->x_H;
    float aH = pow(a, H+1);
    
    float S1;
    float S2;
    if (a != 1 )
    {
        S1 = ( 1 - fabs(a) ) / (1 - fabs(aH)  );
	S2 = ( 1 - fabs(a) ) / ( 1/fabs(aH) - fabs(aH)/aH );
        x->_S1 = S1;				    
    	x->_S2 = S2;				    
    }
}

static void buzz_tilde_ft1(t_buzz_tilde *x, t_float f)
{
    x->x_phase = f;
}


static void buzz_tilde_a(t_buzz_tilde *x, t_float f)
{
    x->x_a = f;
    varsetup(x);
}

static void buzz_tilde_H(t_buzz_tilde *x, t_float f)
{
    x->x_H = f;
    varsetup(x);
}


void buzz_tilde_setup(void)
{
    buzz_tilde_class = class_new(gensym("buzz~"), (t_newmethod)buzz_tilde_new, 0,
    	sizeof(t_buzz_tilde), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(buzz_tilde_class, t_buzz_tilde, x_f);
    class_addmethod(buzz_tilde_class, (t_method)buzz_tilde_dsp, gensym("dsp"), 0);
    class_addmethod(buzz_tilde_class, (t_method)buzz_tilde_ft1, gensym("ft1"), A_FLOAT, 0);
    class_addmethod(buzz_tilde_class, (t_method)buzz_tilde_a, gensym("f_a"), A_FLOAT, 0);
    class_addmethod(buzz_tilde_class, (t_method)buzz_tilde_H, gensym("f_H"), A_FLOAT, 0);
}