Raj and Yoav session in making multiple protein species

Samo Penic

2021-05-10 6487a03362d35589c286d0660eba123a7ab00535

 src/shcomplex.c

@@ -1,8 +1,14 @@
/* vim: set ts=4 sts=4 sw=4 noet : */
#include<math.h>
#include<stdlib.h>
#include<string.h>
#include<gsl/gsl_complex.h>
#include<gsl/gsl_complex_math.h>
#include<gsl/gsl_sf_legendre.h>

#include<gsl/gsl_matrix.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_linalg.h>
#include "general.h"
#include "sh.h"
#include "shcomplex.h"
@@ -64,7 +70,7 @@
    if(sph->co != NULL) free(sph->co);
    if(sph->ulm !=NULL) free(sph->ulm);
    if(sph->ulmComplex !=NULL) free(sph->ulmComplex);

    if(sph->sumUlm2 !=NULL) free(sph->sumUlm2);
        if(sph->Ylmi!=NULL) {
            for(i=0;i<sph->l;i++){
                if(sph->Ylmi[i]!=NULL){
@@ -112,6 +118,27 @@
    return TS_SUCCESS;
}

char *Ulm2Complex2String(ts_vesicle *vesicle){
   ts_int i,j;
   char *strng=(char *)calloc(5000, sizeof(char));
   char tmpstrng[255];
   for(i=0;i<vesicle->sphHarmonics->l;i++){
          for(j=i;j<2*i+1;j++){
         sprintf(tmpstrng,"%e ", gsl_complex_abs2(vesicle->sphHarmonics->ulmComplex[i][j]));
         strcat(strng,tmpstrng);
          }
   }
   //strcat(strng,"\n");

   return strng;
}

ts_bool freeUlm2String(char *strng){
   free(strng);
   return TS_SUCCESS;
}


ts_bool storeUlmComplex2(ts_vesicle *vesicle){

   ts_spharm *sph=vesicle->sphHarmonics;
@@ -125,3 +152,44 @@
   return TS_SUCCESS;
}


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);
    gsl_vector *tau=gsl_vector_alloc(2);
    gsl_vector *b=gsl_vector_alloc(max-min);
    gsl_vector *x=gsl_vector_alloc(2);
    gsl_vector *res=gsl_vector_alloc(max-min);

    //solving (A^T*A)*x=A^T*b
    //fill the data for matrix A and vector b
    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));
            bval=0.0;
            //average for m from 0..l (only positive m's)
            for(j=0;j<=i;j++){
                bval+=vesicle->sphHarmonics->sumUlm2[i][(j+i)];
            }
                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,"b[2]=%e\n",gsl_vector_get(b,1));
    gsl_linalg_QR_decomp(A,tau);
    gsl_linalg_QR_lssolve(A,tau,b,x,res);
//    fprintf(stderr,"kc=%e\n",gsl_vector_get(x,1));
    retval=gsl_vector_get(x,1);
    gsl_matrix_free(A);
    gsl_vector_free(tau);
    gsl_vector_free(b);
    gsl_vector_free(x);
    gsl_vector_free(res);
    
    return retval;
}