Trisurf Monte Carlo simulator
Samo Penic
2014-12-16 6bc95bfe19012a33f185adf1476e623d98d33465
src/shcomplex.c
@@ -9,8 +9,10 @@
#include<gsl/gsl_linalg.h>
#include "general.h"
#include "sh.h"
#include "shreal.h"
#include "shcomplex.h"
#include "string.h"
#include "io.h"
ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l){
    ts_uint j,i;
@@ -51,13 +53,14 @@
    /* Calculate coefficients that will remain constant during all the simulation */ 
   precomputeShCoeff(sph);
   shreal_init(sph,l); //this is for TESTING only!
    return sph;
}
ts_bool complex_sph_free(ts_spharm *sph){
    int i,j;
    if(sph==NULL) return TS_FAIL;
    shreal_free(sph); //this is for TESTING only!
    for(i=0;i<sph->l;i++){
        if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
        if(sph->ulmComplex[i]!=NULL) free(sph->ulmComplex[i]);
@@ -130,9 +133,9 @@
}
ts_double calculateKc(ts_vesicle *vesicle){
    ts_int min=4;
    ts_int max=vesicle->sphHarmonics->l-6;
ts_double calculateKc(ts_vesicle *vesicle, ts_int lmin, ts_int lmax){
    ts_int min=lmin;
    ts_int max=lmax; //vesicle->sphHarmonics->l-3;
    ts_long i,j;
    ts_double retval, bval;
    gsl_matrix *A=gsl_matrix_alloc(max-min,2);
@@ -146,7 +149,7 @@
    for(i=min;i<max;i++){
            gsl_matrix_set(A, i-min,0,(ts_double)((i-1)*(i+2)));
            gsl_matrix_set(A, i-min,1,(ts_double)((i-1)*(i+2)*(i+1)*i));
            fprintf(stderr,"%e %e\n", gsl_matrix_get(A,i-min,0), gsl_matrix_get(A,i-min,1));
//            fprintf(stderr,"%e %e\n", gsl_matrix_get(A,i-min,0), gsl_matrix_get(A,i-min,1));
            bval=0.0;
            //average for m from 0..l (only positive m's)
            for(j=0;j<=i;j++){
@@ -155,7 +158,7 @@
                bval=bval/(ts_double)vesicle->sphHarmonics->N/(ts_double)(i+1);
            gsl_vector_set(b,i-min,1.0/bval);
            fprintf(stderr,"%e\n", 1.0/gsl_vector_get(b,i-min));
//            fprintf(stderr,"%e\n", 1.0/gsl_vector_get(b,i-min));
    }
//    fprintf(stderr,"b[2]=%e\n",gsl_vector_get(b,1));
    gsl_linalg_QR_decomp(A,tau);
@@ -170,3 +173,30 @@
    
    return retval;
}
ts_bool saveAvgUlm2Complex(ts_vesicle *vesicle){
   FILE *fh;
    char filename[10000];
    strcpy(filename, command_line_args.path);
    strcat(filename, "sph2out.dat");
   fh=fopen(filename, "w");
   if(fh==NULL){
      err("Cannot open file %s for writing");
      return TS_FAIL;
   }
   ts_spharm *sph=vesicle->sphHarmonics;
   ts_int i,j;
   fprintf(fh,"l,\tm,\tulm^2avg\n");
   for(i=0;i<sph->l;i++){
          for(j=0;j<2*i+1;j++){
      fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);
          }
    fprintf(fh,"\n");
   }
   fclose(fh);
   return TS_SUCCESS;
}