// SAMPLE C CODE TO ESTIMATE DRAG AND LIFT FROM FLOWSQUARE DATA SET
// COMPILE THIS CODE AND EXECUTE IN THE "DATA" FOLDER.

// This code considers zero-velocity region as object.
// This code assumes a mean convection flow from left to right.
// z-length as unity.
// This code does not account friction drag.

#include <stdio.h>

// Same parameters as in grid.txt
#define nx 512
#define ny 256
#define lx 9.0
#define ly 4.5
// Put your time step
#define fnameu "./ru-000020000.dat"
#define fnamev "./rv-000020000.dat"
#define fnamep  "./p-000020000.dat"

FILE *fpu,*fpv,*fpp;
int i,j;
double ru[ny][nx],rv[ny][nx],p[ny][nx];
double dx,dy,rvelo,px,py;

int main(void){

      dx=lx/(double)(nx-1);
      dy=ly/(double)(ny-1);

      fpu=fopen(fnameu,"rb");
      fread(ru,nx*ny*sizeof(double),1,fpu);
      fpv=fopen(fnamev,"rb");
      fread(rv,nx*ny*sizeof(double),1,fpv);
      fpp=fopen(fnamep,"rb");
      fread(p ,nx*ny*sizeof(double),1,fpp);
      close(fpu);
      close(fpv);
      close(fpp);

      py=0.0;
      px=0.0;
      for(i=0;i<nx;i++){ // Estimating lift
            for(j=0;j<ny;j++){
                  rvelo=ru[j][i]*ru[j][i]+rv[j][i]*rv[j][i];
                  if(rvelo==0.0){
                        py=py+p[j][i];
                        break;
                  }
            }
            for(j=ny-1;j>=0;j--){
                  rvelo=ru[j][i]*ru[j][i]+rv[j][i]*rv[j][i];
                  if(rvelo==0.0){
                        py=py-p[j][i];
                        break;
                  }
            }
      }

      for(j=0;j<ny;j++){ // Estimating drag
            for(i=0;i<nx;i++){
                  rvelo=ru[j][i]*ru[j][i]+rv[j][i]*rv[j][i];
                  if(rvelo==0.0){
                        px=px+p[j][i];
                        break;
                  }
            }
            for(i=nx-1;i>=0;i--){
                  rvelo=ru[j][i]*ru[j][i]+rv[j][i]*rv[j][i];
                  if(rvelo==0.0){
                        px=px-p[j][i];
                        break;
                  }
            }
      }
      px=px*dx;
      py=py*dy;
      printf("Drag: %5e [N]\n",px);
      printf("Lift: %5e [N]\n",py);
}