/*

Principal stress calculator (beta)

Author: Suri Bala
Copyright: Livermore Software
Usage: pstress dynainfile
Compilation: cc -lm pstress.c
Date: Aug 1, 2007

*/


#include<stdio.h>
#include<stdarg.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include<malloc.h>

// definitions

#define LINE_WIDTH 120
#define VAR_WIDTH 20


// prototype declaration

char *getLine(FILE *fp);
float *getPrincipalStress(float sig[]);

// global variables

char tmpbuf[LINE_WIDTH];

int main(int argc,  char *argv[]) 
{

  char *line=NULL;
  FILE *dynain=NULL, *dynaout=NULL;
  int i=0, j=0, eid=0, nplane=0, nthick=0, nhisv=0, avg=1;
  float para=0, avg_sig[6], sig[6], *ps;

  sig[0]=0.0;
  sig[1]=0.0;
  sig[2]=0.0;

  if( argc == 1 ) {
    fprintf(stderr," Usage %s valid_dynain_file\n Exiting ... \n", argv[0]);
    exit(-1);
  }

  dynain = fopen(argv[1], "r");

  if( !dynain ) {
    fprintf(stderr," Could not open dynain file %s\n", argv[1]);
    fprintf(stderr," Exiting ...\n");
    exit(-2);
  }

  dynaout=fopen("dynaout.csv", "w");
  if( !dynaout ) {
    fprintf(stderr," Could not open dynaout file for writing\n");
    fprintf(stderr," Exiting ...\n");
    exit(-3);
  }

  fprintf(dynaout, "eid, s1, s2, s3\n");


  while( !feof(dynain) ) 
  {
    line=getLine(dynain); 
    if( strncmp(line,"*INITIAL_STRESS", 15)==0 ) {
      line=getLine(dynain);
      while( strncmp(line,"*",1)!=0 && !feof(dynain) ) 
      {
        sscanf(line,"%10d%10d%10d%10d", &eid, &nplane, &nthick, &nhisv);
	avg=1.0;
        for(i=0; i<nthick; i++ ) 
	{
	  sig[0]=sig[1]=sig[2]=sig[3]=sig[4]=sig[5]=0;
	  avg_sig[0]=avg_sig[1]=avg_sig[2]=avg_sig[3]=avg_sig[4]=avg_sig[5]=0;
          for(j=0; j<nplane; j++) 
	  {
            line=getLine(dynain);
            sscanf(line,"%10f%10f%10f%10f%10f%10f%10f", &para, &sig[0], &sig[1], &sig[2], &sig[3], &sig[4], &sig[5]);
	    avg_sig[0] += sig[0];
	    avg_sig[1] += sig[1];
	    avg_sig[2] += sig[2];
	    avg_sig[3] += sig[3];
	    avg_sig[4] += sig[4];
	    avg_sig[5] += sig[6];
	    if( sig[0]>0 && sig[1]>0 && sig[2] && sig[3]>0 && sig[4]>0 && sig[5]>0 ) ++avg;
   	  }
  	  if( i==0 || i==nthick-1) {
	    avg_sig[0] = avg_sig[0]/avg;
	    avg_sig[1] = avg_sig[1]/avg;
	    avg_sig[2] = avg_sig[2]/avg;
	    avg_sig[3] = avg_sig[3]/avg;
	    avg_sig[4] = avg_sig[4]/avg;
	    avg_sig[5] = avg_sig[5]/avg;
  	    ps=getPrincipalStress(avg_sig);
  	    fprintf(dynaout,"%d, %10f, %10f, %10f\n", eid, ps[0], ps[1], ps[2]);
  	  }
        }
        line=getLine(dynain);
      }
    }
  }

  fclose(dynain);

}


char *getLine(FILE *fp) {
 int i=0, len=0;
 fgets(tmpbuf, LINE_WIDTH, fp);
 while( strncmp(tmpbuf,"$",1)==0 ) {
  fgets(tmpbuf, LINE_WIDTH, fp);
 }
 len=strlen(tmpbuf);
 tmpbuf[len-1]='\0';
 for(i=0; i<len; i++){
   tmpbuf[i] = toupper(tmpbuf[i]);
 }
 return tmpbuf;
}

float *getPrincipalStress(float sig[]) {
 
    float iv[3], *ps=NULL;
    float f=0, g=0, h=0, r=0, s=0, t=0, u=0, l=0, m=0, n=0, p=0, a=0, b=0, c=0, d=0, i=0,j=0,k=0; 


    ps=(float *)malloc(3*sizeof(float));

    iv[0] = sig[0]+sig[1]+sig[2];
    iv[1] = sig[0]*sig[1]+sig[1]*sig[2]+sig[2]*sig[1]-pow(sig[3],2)-pow(sig[4],2)-pow(sig[5],2);
    iv[2] = sig[0]*sig[1]*sig[2]+2*sig[3]*sig[4]*sig[5]-sig[0]*pow(sig[4],2)-sig[1]*pow(sig[5],2)-sig[2]*pow(sig[3],2);

    a=1.0;
    b=-iv[0];
    c=iv[1];
    d=-iv[2];

    f = ( (3.0*c/a) - (pow(b,2)/pow(a,2.0)) ) / 3.0;
    g = ( (2.0*pow(b,3.0)/pow(a,3.0)) - (9.0*b*c/pow(a,2.0) ) + (27.0*d/a) ) / 27.0;
    h = ( (pow(g,2.0)/4.0) + (pow(f,3.0)/27.0) );
    i = sqrt((pow(g,2.0)/4.0)-h);
    j=pow(i, 1.0/3.0);
    k = acos(-(g/(2.0*i) ));

    if( h < 0 ) {
       // all three roots are real 
       l=j*-1.0;
       m=cos(k/3.0);
       n=sqrt(3.0)*sin(k/3.0);
       p=(b/3.0*a)*-1.0;

       ps[0]=2.0*j*cos(k/3.0)-(b/(3.0*a));
       ps[1]=l*(m+n)+p;
       ps[2]=l*(m-n)+p;
    } else if(h==0 && g==0 && f==0) {
      ps[0] = ps[1] = ps[2] = pow(d/a, 1/3)*-1;
    } else if(h>0 ) {
       r=-(g/2.0)+sqrt(h);
       s=pow(r, 1/3);
       t=-(g/2)-sqrt(h);
       u=pow(t,1/3);

       ps[0] = (s+u)-(b/3*a);
       ps[1] = -((s+u)/2) - (b/3*a) + i*(s-u)*sqrt(3)*0.5;
       ps[2] = -((s+u)/2) - (b/3*a) - i*(s-u)*sqrt(3)*0.5;
    }

    return ps;

}