#include<math.h>

void lpE( double *E, double *theta, double *t2s, double *s2, int *n, int *m, 
          double *sm)
{
   int i, j; 
   double the, se, mj, vj ; 
   for(j = 0 ; j < *m ; j++ ) {
       vj=1.0/( (*n)*1.0/s2[j] +1.0/(*t2s*s2[j]) ); 
       se=0.0;
       for(i=0 ; i< *n ; i++ ) { se=se + E[ i*(*m) + j ] ; }
       mj=vj*se/s2[j] ; 
       the= .5*(log(vj) - log(*t2s) - log(s2[j]) +mj*mj/vj );
       *sm=*sm+ the +log( exp(-the)+exp(*theta ))- log(1.0+exp( *theta)) ;
                              }
}