gaussian_kernel.c

#include "gaussian_kernel.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
struct gaussian_kernel
{
    double sigma;
    double amp;
    double cutoff;
};

void* gaussian_create(char* name)
{
    void* src = malloc(sizeof(gaussian_kernel));
    memset(src, 0, sizeof(gaussian_kernel));
    FILE* infile = fopen(name, "r");
    if(infile == NULL)
    {
        perror("Error in opening file:");
        free(src);
        return NULL;
    }
    char buffer[512];
    while(fgets(buffer, 512, infile) != NULL)
    {
         char* tmp = strtok(buffer, " ");
         if(strcmp(tmp, "Sigma") == 0)
         {
             tmp = strtok(NULL, " ");
             if(tmp != NULL)
                 ((gaussian_kernel*)src)->sigma = atof(tmp);
             else
             {
                 printf("Error formating\n");
                 free(src);
                 return NULL;
             }
         }
         else if(strcmp(tmp, "Amp") == 0)
         {
             tmp = strtok(NULL, " ");
             if(tmp != NULL)
                 ((gaussian_kernel*)src)->amp = atof(tmp);
             else
             {
                 printf("Error formating\n");
                 free(src);
                 return NULL;
             }
         }
         else if(strcmp(tmp, "Cutoff") == 0)
         {
             tmp = strtok(NULL, " ");
             if(tmp != NULL)
                 ((gaussian_kernel*)src)->cutoff = atof(tmp);
             else
             {
                 printf("Error formating\n");
                 free(src);
                 return NULL;
             }
         }
         else
         {
             printf("Unrecognized pattern\n");
             free(src);
             return NULL;
         }
    }
    fclose(infile);
    return src;
}

double gaussian_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
{
    gaussian_kernel* kernel = (gaussian_kernel*)src;
    double pos[3];
    for(size_t dim = 0; dim < 3; ++dim)
        pos[dim] = basis[dim][0]*i + basis[dim][1]*j + basis[dim][2]*k + origin[dim];
    double d = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2];    
    if(d > (kernel->cutoff)*(kernel->cutoff))
        return 0;
    else
        return kernel->amp*exp(-d/((kernel->sigma)*(kernel->sigma)));
}

void gaussian_destroy(void* src)
{
    free(src);
}