ComplianceToModuli.c

/*

---------------------------------------------------------------------------

  C program to convert comliance data to storage and loss moduli, using the

  method described in Phys. Rev. E 80, 012501 (2009).

R M L Evans 6/10/2009

---------------------------------------------------------------------------

*/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

void main(void)

{

  FILE *fp;

  double t[20000],J[20000],Gp[1000],Gpp[1000];

  double J0,eta,ReSum,ImSum,x,y,w,wmin,wmax,s,smin,smax,ds;

  int N,j,k,nw;

  char filename1[100],filename2[100];

  char answer;

  printf("Enter name of ascii file containing (time,compliance) data:\n");

  scanf("%s",filename1);

  printf("Enter name of output file:\n");

  scanf("%s",filename2);

  if (fp = fopen(filename1,"r"))

  {

    N=1;

    while (fscanf(fp,"%lf %lf", &t[N], &J[N]) > 0) N++;   //Read data until end of set

    fclose(fp);

  }

  else exit(1);

  N--;

  printf("Zero-time extrapolation of compliance: J0=");

  scanf("%lf",&J0);

  printf("Long-time viscosity: eta=");

  scanf("%lf",&eta);

  printf("Number of frequencies: ");

  scanf("%d",&nw);

  while ((answer != 'Y')&&(answer != 'y')&&(answer != 'N')&&(answer !='n'))

  {

    printf("\nAngular frequency range from 1/t_Max to 1/t_1 ? (Y/N): ");

    scanf("%c",&answer);

  }

  if ((answer == 'Y')||(answer == 'y'))

  {

    wmin = 1./t[N];

    wmax = 1./t[1];

  }

  else

  {

    printf("\nMinimum angular frequency: omega_Min=");

    scanf("%lf",&wmin);

    printf("\nMaximum angular frequency: omega_Max=");

    scanf("%lf",&wmax);

  }

  smin = log10(wmin);

  smax = log10(wmax);

  ds = (smax-smin)/(nw-1);

  for (j=0; j<nw; j++)

  {

    s = smin + (j*ds);

    w = pow(10.0,s);

    ReSum = 0.;

    ImSum = 0.;

    for (k=2; k<=N; k++)

    {

      ReSum += (cos(w*t[k-1])-cos(w*t[k])) * (J[k]-J[k-1])/(t[k]-t[k-1]);

      ImSum += (-sin(w*t[k-1])+sin(w*t[k])) * (J[k]-J[k-1])/(t[k]-t[k-1]);

    }

    x = ((1.-cos(w*t[1]))*(J[1]-J0)/t[1]) + (cos(w*t[N])/eta) + ReSum;

    y = (w*J0) + (sin(w*t[1])*(J[1]-J0)/t[1]) - (sin(w*t[N])/eta) + ImSum;

    Gp[j] = w*y/((x*x)+(y*y));

    Gpp[j] = w*x/((x*x)+(y*y));

  }

  if (fp = fopen(filename2,"w"))

  {

    for (j=0; j<nw; j++) fprintf(fp,"%e  %e  %e\n",pow(10.0,smin+(j*ds)),Gp[j],Gpp[j]);

    fclose(fp);

  }

  printf("\nOutput written to %s\n",filename2);

  printf("\nPress any key to finish.\n");

  getch();

  return;

}