#ifndef _BEAM_H_
#define _BEAM_H_


#define NSTRN 6
#define NDIM 3
#define MEMB_CONST 7
#define NDOF_ND 12

#define CHKINP(r) if (r==-1){ printf("Error reading input file (line %d)\n",__LINE__); exit(EXIT_FAILURE);}
#define CHKFRMT(r,val) if (r!=val){ printf("Error formatting input file (line: %d)\n",__LINE__); exit(EXIT_FAILURE);}

#define ALLOCATE1(ptr,type,n1) ptr = (type*)malloc((n1)*sizeof(type)) 
#define ALLOCATE2(ptr,type,n1,n2) ptr = (type **)malloc((n1) * sizeof(type *)); \
                                  for (i_=0; i_<(n1); i_++){\
                                    ptr[i_] = (type *)malloc((n2) * sizeof(type));}
#define ALLOCATE3(ptr,type,n1,n2,n3) ptr = (type ***)malloc(((n1)) * sizeof(type **)); \
                                  for (i_=0; i_<(n1); i_++){\
                                    ptr[i_] = (type **)malloc((n2) * sizeof(type *));\
                                    for (j_=0; j_<(n2); j_++){\
                                        ptr[i_][j_] = (type *)malloc((n3) * sizeof(type));}}

#define DEALLOCATE1(ptr) free(ptr);

#define DEALLOCATE2(ptr,n1) for (i_=0; i_<(n1); i_++){\
                               free(ptr[i_]);}\
                            free(ptr);

#define DEALLOCATE3(ptr,n1,n2) for (i_=0; i_<(n1); i_++){\
                        for (j_=0; j_<(n2); j_++){\
                           free(ptr[i_][j_]);}\
                           free(ptr[i_]);}\
                           free(ptr);


#define READL(line,len,file) read = getline(&line,&len,file); CHKINP(read);\
                                    if(line==NULL){printf("Reached end of the file\n");return;}\
                                    while (line[0]=='\n' || line[0]=='\r' || line[0]=='#'){\
                                        read = getline(&line,&len,file); CHKINP(read);\
                                        if(line==NULL){printf("Reached end of the file\n");return;}\
                                    }

#define READ1(file,val,tag) READL(line,len,file)\
                        sscanf(line,#tag,&(val))

#define READ3(file,vec,tag) READL(line,len,file)\
                        sscanf(line,#tag #tag #tag,&vec[0],&vec[1],&vec[2])

                        


struct PrescriInf{
	 int   id;                               // where it is applied, could be a node number or member number
	 int   dof[NSTRN];                       // maximum 6 degrees of freedom can be prescribed, 	                         
	                                         // for distributed loads, it is used to denote the distribution function no
     double value[NSTRN];                    // the magnitude of the prescribed values
	 int    time_fun_no[NSTRN];              // which time function is used
	 int follower[NSTRN];                    // whether the prescribed quantity is a follower or not: 1 is a follower; 0 is not
	 double value_current[NSTRN];            // indicate the current functional value updated by time steps, calculated internally
};

struct TimeFunction{
	 int fun_type;                           //unction type: 0, user defined; 1, harmonic
	 double ts, te;                          // starting and ending time
	 int entries;                            // number of entries
	 double *time_val;                       // the ith time, in increasing order; the amplitude of the harmonic
	 double *fun_val ;                       // the ith functional value; the period of the harmonic
	 double *phase_val;                      // the phase value of the harmonic
};





typedef struct
{
    // Declare all the variables required for the analysis
    int                nkp;
    int                nelem;             // number of elements
    int                ndof_el;           // number of unknowns for each element
    int                nmemb;             // number of beam members
    int                ncond_pt;          // number of point conditions
    int                nmate;             // number of different cross-sections
    int                nframe;            // number of different cross-sectional frames
    int                ndistrfun;         // number of distribution functions
    int                ncurv;             // number of sets of initial curvature/twists
    double             **coord;           // coordinates for key points
    int                **member;          // member desriptions
    struct PrescriInf  *pt_condition;     // prescribed concentrated infromation at points
    double             ***material;       // cross-sectional properties

    int                ntimefun;          // number of time functions 
    int                niter;             // maximum number of iterations
    int                nstep;             // number of time steps
    double             ***sol_pt;         // solutions for points
    double             ***sol_mb;         // solutions for members
    char               error[32];            // error message

    int                ncond_mb;         // number of prescribed distribution               			                		
    double             ***frame;          // the direction cosine matrix for each frame. The default frame is the global frame 
    struct PrescriInf  *mb_condition;     // distributed load information
    double             **distr_fun;       // distributed functions 
    double             **curvature;       // curvature information
    double             omega_a0[NDIM];    // angular velocity of frame a
    double             v_root_a0[NDIM];   // linear velocity of starting point of the first member
    int                omega_a_tf[NDIM];  // time function for angular velocity of frame a
    int                v_root_a_tf[NDIM]; // time function for linear velocity of starting point of the first member

    double               simu_time[2];    // start and end time
    struct TimeFunction  *time_function;  // time functions
    double               **init_cond;
    int                  analysis_flag;
    int                  nev;             // upon return, nev is the number of converged eigen values might be the original nev given by the user in the inputs or nev+1 due to complex conjugate pair 

    // Personal add
    double              L;                // Beam Length
    int                 initFromFile;
} Beam;



void initBeam(Beam *bm);
void beam_setLoads(Beam *bm, double *loads, int load_no);
void beam_writeSolToFile(Beam *bm);
void freeBeam(Beam *bm);



#endif