/********************************************************************** 
 * Filename: qjairflow.ch                                             *
 * Description: This program is coverted from a fortran program       * 
 *    qjairflow.f used in EAE127.  To run this program, use computers *
 *    in CAE lab, or download ChSciTE IDE from                        *
 *    http://chscite.sourceforge.net and C/C++ interpreter Ch         *
 *    Professional or Student Edition from                            *
 *     http://www.softintegration.com/download                        *
 ----------------------------------------------------------------------
 History:                                                           
  *) 6/2/07,  Converted by 
     Yu-Cheng Chou cycchou@ucdavis.edu                            
**********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <numeric.h>

#define ix 41

void moran(int iprv, int ixx, int &itrans, int &isep, double Re, 
           array double x[1:ix], array double y[1:ix], array double ve[1:ix], 
           array double xx[1:ix], array double vgrad[1:ix], 
           array double theta[1:ix], array double cfr[1:ix]);
void thwats(double alam, double &h, double &el);
void runge2(int ixx, int i, double Re, double dx, array double yy[1:2],
            array double ve[1:ix], array double vgrad[1:ix]);
void derivs(int ixx, int i, double Re, array double yp[1:2],
            array double yt[1:2], array double ve[1:ix],
            array double vgrad[1:ix]); 
double h1ofh(double h);
double hofh1(double h1);
double cfturb(double rtheta, double h);

int main() {
    array double x[1:ix], axi[1:ix], zu[1:ix], zo[1:ix];
    array double f[1:ix], fp[1:ix];
    array double e[1:ix], ep[1:ix];
    array double cpu[1:ix], cpo[1:ix];
    array double vu[1:ix], vo[1:ix];
    array double cpusd[1:ix], cposd[1:ix];
    array double ve[1:ix], xx[1:ix], vgrad[1:ix], theta[1:ix], cfr[1:ix];
    double complex un, im, zl, dzds, dfdzu, dfdzo;
    int itrans, isep;
    FILE *fp1, *fp2, *fp3, *fp4, *fp5, *fp6, *fp7, *fp8, *fp9;

    fp1 = fopen("qjairflow.dat", "r");
    fp2 = fopen("qjairflow.xzu", "w");
    fp3 = fopen("qjairflow.xzo", "w");
    fp4 = fopen("qjairflow.cpu", "w");
    fp5 = fopen("qjairflow.cpo", "w");
    fp6 = fopen("qjairflow.cfu", "w");
    fp7 = fopen("qjairflow.cfo", "w");
    fp8 = fopen("qjairflow.cpusd", "w");
    fp9 = fopen("qjairflow.cposd", "w");

    double eps = 1.e-7;
    double pi = 2.*asin(1.);
    double degrad = pi/180.;
    double dtet = pi/(ix-1);
    un = complex(1., 0.);
    im = complex(0., 1.);

    double Uref = 1.;
    double dm, em, alphd, delt, epsi, alpha, sina, cosa, teta;

    fscanf(fp1, "%lf%*s%*s%*s%*s%*s", &dm);
    fscanf(fp1, "%lf%*s%*s%*s%*s", &em);
    fscanf(fp1, "%lf", &alphd);
    fclose(fp1);

    //printf("dm = %lf\n", dm);
    //printf("em = %lf\n", em);
    //printf("alphd = %lf\n", alphd);

    delt = 2.*dm;
    epsi = 4./(3.*sqrt(3.))*em+1.0e-3;
    alpha = alphd*degrad;
    sina = sin(alpha);
    cosa = cos(alpha);
    teta = -dtet;

    //printf("delt = %lf\n", delt);
    //printf("epsi = %lf\n", epsi);
    //printf("alpha = %lf\n", alpha);
    //printf("sina = %lf\n", sina);
    //printf("cosa = %lf\n", cosa);
    //printf("teta = %lf\n", teta);
    
    int i, ic;
    double sint, cost, xi, zui, zoi, cpui, cpoi, anorm;
    double complex z1, tmp;

    for(i=1; i<=ix; i++) {
        ic = ix+1-i;

        //if(ic == 1) {
        //    printf("teta = %f\n", teta);
        //    printf("dtet = %f\n\n", dtet);
        //}

        teta = teta+dtet;
        sint = sin(teta);
        cost = cos(teta);
        xi =.5*(1.+cost);
        x[ic] = xi;
        zui = .5*(epsi*(1.-cost)*sint+delt*sint*sint);
        zoi = .5*(-epsi*(1.-cost)*sint+delt*sint*sint);
        zu[ic] = zui;
        zo[ic] = zoi;
        z1 = (1.+epsi*(1.-cost)+delt*sint)*exp(im*teta);
        
        //if(ic == 1) {
        //    printf("epsi= %f\n", epsi);
        //    printf("cost= %f\n", cost);
        //    printf("delt= %f\n", delt);
        //    printf("sint= %f\n", sint);
        //    printf("im= %f\n", im);
        //    printf("teta= %f\n", teta);
        //    printf("z1= %f\n", z1);
        //}

        //printf("%f\n", z1);

        tmp = z1+un*epsi-im*delt;
        tmp *= tmp;
        dfdzu = (exp(-im*alpha)+2.*(im*sina+un*(epsi*cosa-delt*sina))
                 /(z1+1.))*(z1*z1/tmp);

        z1 = (1.+epsi*(1.-cost)-delt*sint)*exp(-im*teta);

        //if(ic == 1) {
        //    printf("\n");
        //    printf("epsi= %f\n", epsi);
        //    printf("cost= %f\n", cost);
        //    printf("delt= %f\n", delt);
        //    printf("sint= %f\n", sint);
        //    printf("im= %f\n", im);
        //    printf("teta= %f\n", teta);
        //    printf("z1= %f\n", z1);
        //}

        //printf("%f\n", z1);

        tmp = z1+un*epsi-im*delt;
        tmp *= tmp;
        dfdzo = (exp(-im*alpha)+2.*(im*sina+un*(epsi*cosa-delt*sina))
                 /(z1+1.))*(z1*z1/tmp);

        //printf("%f    %f\n", dfdzu, dfdzo);

        cpui = -1.+abs(dfdzu)*abs(dfdzu);
        cpoi = -1.+abs(dfdzo)*abs(dfdzo);
        cpu[ic] = cpui;
        cpo[ic] = cpoi;
        dzds = un*(2.*sint)-im*(2.*epsi*(cost-cost*cost+sint*sint)
                                +4.*delt*sint*cost);
        anorm = abs(dzds);

        if(anorm < eps) {
            dzds = (un-im*2.*delt)/sqrt(1.+4.*delt*delt);
        }
        else {
            dzds = dzds/(anorm+eps);
        }

        vu[ic] = real(dfdzu*dzds);

        //if(ic == 1) {
        //    printf("dfdzu= %f  dzds= %f  vu[1]= %f\n", dfdzu, dzds, vu[1]);
        //}

        dzds = un*(2.*sint)+im*(2.*epsi*(cost-cost*cost+sint*sint)
                                -4.*delt*sint*cost);
        anorm=abs(dzds);

        if(anorm < eps) {
            dzds = (un-im*2.*delt)/sqrt(1.+4.*delt*delt);
        }
        else {
            dzds = dzds/(anorm+eps);
        }

        vo[ic] = real(dfdzo*dzds);

        //if(ic == 1) {
        //    printf("dfdzo= %f  dzds= %f  vo[1]= %f\n", dfdzo, dzds, vo[1]);
        //}

        cpusd[ic] = 8.*em*(1.-4.*xi/3.)/sqrt(3.)
                    +2.*alpha*sqrt((1.-xi)/(xi+eps))
                    +16.*dm*sqrt(xi*(1.-xi));
        cposd[ic] = 8.*em*(1.-4.*xi/3.)/sqrt(3.)
                    -2.*alpha*sqrt((1.-xi)/(xi+eps))
                    -16.*dm*sqrt(xi*(1.-xi));

        //printf("cpusd= %f    cposd= %f\n", cpusd[ic], cposd[ic]);
    }

    //printf("vu[1] = %f    vo[1] = %f\n", vu[1], vo[1]);

    int ipr;
    double axii, fim, eim, axiim, teti;
 
    printf("*************************\n");
    printf("dimensionless parameters:\n");
    printf("                           X=C*x\n");
    printf("                           Z=C*z\n");
    printf("                      P-Pinf=0.5*Rho*V**2*Cp\n");
    printf("                           L=0.5*Rho*V**2*C*Cl\n");
    printf("                           D=0.5*Rho*V**2*C*Cd\n");
    printf("                        CM,o=0.5*Rho*V**2*C**2*Cm,o\n");
    printf("\n"); 
    printf("*************\n");
    printf("profile data:\n");
    printf("             relative camber= %f\n", dm);
    printf("          relative thickness= %f\n", em);
    printf("             angle of attack= %f (rd) = %f (deg)\n", alpha, alphd);
    printf("do you want to print the mesh, geometry and flow? y/n=1/0\n");
    scanf("%d", &ipr);

    if(ipr == 1) {
        printf("\n");    
        printf("  i=     x(i)=     f(i)=    fp(i)=     "\
               "e(i)=    ep(i)=    vu(i)=    vo(i)=\n");
    }

    axii = .5*(1.-cos(.5*dtet));
    fim = -4.*dm*axii*(1.-axii);
    eim = -.5*em*(1.+cos(.5*dtet))*sin(.5*dtet);
    axiim = -axii;

    for(i=1; i<=ix; i++) {
        teti = (i-1)*dtet;
        xi = .5*(1.-cos(teti));
        axii = .5*(1.-cos(teti+.5*dtet));
        x[i] = xi;
        axi[i] = axii;
        f[i] = 4.*dm*axii*(1.-axii);
        fp[i] = (f[i]-fim)/(axii-axiim);
        e[i] = .5*em*(1.+cos(teti+.5*dtet))*sin(teti+.5*dtet);
        ep[i] = (e[i]-eim)/(axii-axiim);
        eim = e[i];
        e[i] = .5*em*(1.+cos(teti))*sin(teti);
        fim = f[i];
        f[i] = 4.*dm*xi*(1.-xi);
        axiim = axii;

        if(i == (ix-1)) {
            fim = -fim;
            eim = -3.*eim;
            axiim = 2.-axii;
        }

        if(ipr == 1) {
            printf(" %3d  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %9.4f\n", 
                   i, x[i], f[i], fp[i], e[i], ep[i], vu[i], vo[i]);
        }
        
        fprintf(fp2, "%f  %f\n", x[i], zu[i]);
        fprintf(fp3, "%f  %f\n", x[i], zo[i]);
        fprintf(fp4, "%f  %f\n", x[i], cpu[i]);
        fprintf(fp5, "%f  %f\n", x[i], cpo[i]);
        fprintf(fp8, "%f  %f\n", x[i], cpusd[i]);
        fprintf(fp9, "%f  %f\n", x[i], cposd[i]);
    }

    fclose(fp2);
    fclose(fp3);
    fclose(fp4);
    fclose(fp5);
    fclose(fp8);
    fclose(fp9);
    
    double Cl, Cd, Cmo, Clsd, Cdsd, Cmosd, Re, Cdt;
    int iprv;

    printf("\n");
    printf("****************************************\n");
    printf("exact inviscid aerodynamic coefficients:\n");

    Cl = 2.*pi*((1.+epsi)*sina+delt*cosa);
    Cd = 0.;
    printf("Cl= %f    Cd= %f\n\n", Cl, Cd);

    Cmo = -.5*pi*((1.+epsi)*cosa*sina+2.*delt);
    printf("Cm,o= %f\n\n", Cmo); 

    printf("*******************************************\n");
    printf("small disturbance aerodynamic coefficients:\n");

    Clsd = 2.*pi*(alpha+delt);
    Cdsd = 0.;
    printf("Clsd= %f   Cdsd= %f\n\n", Clsd, Cdsd); 

    Cmosd = -.5*pi*(alpha+2.*delt);
    printf("Cm,osd= %f\n\n", Cmosd); 

    printf("          ************************\n");
    printf("**********BOUDARY LAYER CORRECTION************\n");
    printf("          ************************\n\n");

    while(1) {
        printf("Reynolds Number Re=? (Re<1.E-7 or Re>1.E7 to exit)\n");
        scanf("%lf", &Re);

        if((Re<=eps) || (Re>=1./eps)) {
            printf("\n");
            printf("******input files:\n");
            printf("qjairflow.dat                    :dm,em,alpha\n");
            printf("******output files:\n");
            printf("qjairflow.xzu                    :x, zu(x)\n");
            printf("qjairflow.xzo                    :x, zo(x)\n");
            printf("qjairflow.cpu                    :x, Cpu(x)\n");
            printf("qjairflow.cpo                    :x, Cpo(x)\n");
            printf("qjairflow.cfu                    :x, cfu(x)\n");
            printf("qjairflow.cfo                    :x, cfo(x)\n");
            printf("qjairflow.cpusd                  :x, Cpusd(x)\n");
            printf("qjairflow.cposd                  :x, Cposd(x)\n");

            fclose(fp6);
            fclose(fp7);

            return 0;
        }

        printf("\ndo you want to print the viscous flow? y/n=1/0\n");
        scanf("%d", &iprv);
 
        if(iprv == 1) {
            printf("\nupper surface\n");
        }

        for(i=1; i<=ix; i++) {
            ve[i] = vu[i];
            f[i] = f[i]+e[i];
        }

        moran(iprv, ix, itrans, isep, Re, x, f, ve, xx, vgrad, theta, cfr);

        Cd = 0.;

        for(i=2; i<=ix; i++) {
            Cd = Cd+.5*(cfr[i-1]+cfr[i])*(x[i]-x[i-1]);
            if(i > isep) {
                cfr[i] = 0.;
            }
            fprintf(fp6, "%f  %f\n", x[i], cfr[i]);
        }

        printf("upper surface itrans= %d  isep= %d  Cd= %f\n", itrans, isep, Cd);
        Cdt = Cd;
 
        if(iprv == 1) {
            printf("\n\nlower surface\n");
        }
    
        for(i=1; i<=ix; i++) {
            ve[i] = vo[i];
            f[i] = f[i]-2.*e[i];
        }
 
        moran(iprv, ix, itrans, isep, Re, x, f, ve, xx, vgrad, theta, cfr);

        Cd = 0.;

        for(i=2; i<=ix; i++) {
            Cd = Cd+.5*(cfr[i-1]+cfr[i])*(x[i]-x[i-1]);
            if(i > isep) {
                cfr[i] = 0.;
            }
            fprintf(fp7, "%f  %f\n", x[i], cfr[i]);
        }

        printf("lower surface itrans= %d  isep= %d  Cd= %f\n", itrans, isep, Cd);
 
        Cdt = Cdt+Cd;

        printf("\nCd total= %f\n", Cdt);
    }

    fclose(fp6); 
    fclose(fp7); 
 
    return 0;
}

void moran(int iprv, int ixx, int &itrans, int &isep, double Re, 
           array double x[1:ix], array double y[1:ix], array double ve[1:ix], 
           array double xx[1:ix], array double vgrad[1:ix], 
           array double theta[1:ix], array double cfr[1:ix]) {
    int i;
    double dx, dy, v0, v1, x1, v2, x2, v3, x3, fact, h, cf, cfe, dv;
    double alam, el, dth2ve6, Rex, Ret, Retmax, rtheta;
    array double yy[1:2];

    itrans = ixx;
    isep = ixx;
    double pi = 3.1415926535;
    double eps = 1.e-12;
    xx[1] = 0.;

    for(i=2; i<=ixx; i++) {
        dx = x[i]-x[i-1];
        dy = y[i]-y[i-1];
        xx[i] = xx[i-1]+sqrt(dx*dx+dy*dy);
    }

    v1 = ve[3];
    x1 = xx[3];
    v2 = ve[1];
    x2 = xx[1];

    for(i=1; i<=ixx; i++) {
        v3 = v1;
        x3 = x1;
        v1 = v2;
        x1 = x2;

        if(i < ixx) {
            v2 = ve[i+1];
            x2 = xx[i+1];
        }
        else {
            v2 = ve[ixx-1];
            x2 = xx[ixx-1];
            v3 = ve[ixx-2];
            x3 = xx[ixx-2];
        }
        
        fact = (x3-x1)/(x2-x1);
        vgrad[i] = ((v2-v1)*fact-(v3-v1)/fact)/(x3-x2);
    }

    if(iprv == 1) {
        printf("\n       x         y           ve           vdot         theta            h            "\
               "cf       cf-flat plate\n\n");
    }

    i = 1;

    while(1) {
        v0 = ve[i];
        while(ve[i] < 0.) {
            theta[i] = 0.;
            h = 0.;
            cf = 0.;
            cfr[i] = cf;
            cfe = 0.;
            if(iprv == 1) {
                printf("l  %8.4f  %8.4f  %12.4e  %12.4e  %12.4e  %12.4e  %12.4e  %12.4e\n", 
                        x[i], y[i], ve[i], vgrad[i], theta[i], h, cf, cfe);
            }
            i = i+1;
            v0 = ve[i];
        }

        dv = eps*vgrad[i];
        if(dv < pow(eps, 2)) {
            dv = pow(eps, 2);
        }

        theta[i] = sqrt(eps*(6.*pow(v0, 5)+15.*pow(v0, 4)*dv+20.*pow(v0, 3)*
                    pow(dv, 2)+15.*pow(v0, 2)*pow(dv, 3)+6.*v0*pow(dv, 4)+
                    pow(dv, 5))/(6.*pow((v0+dv), 6)));
        do {
            alam = pow(theta[i], 2)*vgrad[i]*Re;
            if(alam > .1) {
                alam =.1;
            }
            if(alam < -.0841) {
                alam =-.0841;
            }

            /***************************************************** 
             * The 'if statement' below is unnecessary. However, * 
             * since it appears in the original Fortran code, I  * 
             * still did the conversion for this part.           * 
             *****************************************************/ 
            if(alam < -.0842) {
                isep = i;
                if(iprv == 1) {
                    printf("laminar separation\n");
                }
                return;
            }

            thwats(alam, h, el);

            if(ve[i] < eps) {
                cf = 2.*el/Re/theta[i]/eps;
            }
            else {
                cf = 2.*el/Re/theta[i]/ve[i];
            }

            cfr[i] = cf;
            cfe = .664/sqrt(Re*(x[i]+eps-x[1])/(x[ixx]-x[1]));

            if(iprv == 1) {
                printf("l  %8.4f  %8.4f  %12.4e  %12.4e  %12.4e  %12.4e  %12.4e  %12.4e\n", 
                        x[i], y[i], ve[i], vgrad[i], theta[i], h, cf, cfe);
            }

            i = i+1;
            if(i > ixx) {
                return;
            }

            if(ve[i] < 0.) {
                break;
            }
            dth2ve6 = .225*(pow(ve[i], 5)+pow(ve[i-1], 5))*(xx[i]-xx[i-1])/Re;
            theta[i] = sqrt((pow(theta[i-1], 2)*pow(ve[i-1], 6)+dth2ve6)/pow(ve[i], 6));
            Rex = Re*xx[i]*ve[i];
            Ret = Re*theta[i]*ve[i];
            Retmax = 1.174*(1.+22400./Rex)*pow(Rex, .46);

        } while((Ret<Retmax) || (i <= 1));

        if(ve[i] >= 0.) {
            break;
        }
    }
    
    itrans = i;
    h = 1.47;
    i = itrans;
    yy[2] = h1ofh(h);
    yy[1] = theta[i-1];

    do {
        dx = xx[i]-xx[i-1];

        runge2(ixx, i, Re, dx, yy, ve, vgrad); 

        theta[i] = yy[1];
        h = hofh1(yy[2]);
        rtheta = Re*ve[i]*theta[i];
        cf = cfturb(rtheta, h);
        cfr[i] = cf;
        cfe = .0592/pow(Re*(x[i]+eps-x[1])/(x[ixx]-x[1]), .2);

        if(iprv == 1) {
            printf("t  %8.4f  %8.4f  %12.4e  %12.4e  %12.4e  %12.4e  %12.4e  %12.4e\n", 
                    x[i], y[i], ve[i], vgrad[i], theta[i], h, cf, cfe);
        }

        if(h > 2.4) {
            isep = i;
            if(iprv == 1) {
                printf("turbulent separation\n");
            }
            return;
        }

        i = i+1;

    } while(i <= ixx);
}

void thwats(double alam, double &h, double &el) { 
    if(alam < 0.) {
        el = .22+1.402*alam+.018*alam/(.107+alam);
        h = 2.088+.0731/(.14+alam);
    }
    else {
        el = .22+alam*(1.57-1.8*alam);
        h = 2.61-alam*(3.75-5.24*alam);
    }
} 

void runge2(int ixx, int i, double Re, double dx, array double yy[1:2], 
            array double ve[1:ix], array double vgrad[1:ix]) {
    int j;
    array double ys[1:2], yp[1:2], yt[1:2];

    for(j=1; j<=2; j++) {  
        yt[j] = yy[j];
    }

    derivs(ixx, i-1, Re, yp, yt, ve, vgrad);

    for(j=1; j<=2; j++) {
        yt[j] = yy[j]+dx*yp[j];
        ys[j] = yy[j]+.5*dx*yp[j];
    }

    derivs(ixx, i, Re, yp, yt, ve, vgrad);

    for(j=1; j<=2; j++) {
        yy[j] = ys[j]+.5*dx*yp[j];
    }
}

void derivs(int ixx, int i, double Re, array double yp[1:2], 
            array double yt[1:2], array double ve[1:ix], 
            array double vgrad[1:ix]) {
    double h1, h, rtheta;

    h1 = h1ofh(hofh1(yt[2]));
    h = hofh1(h1);
    rtheta = Re*ve[i]*yt[1];
    yp[1] = -(h+2.)*yt[1]*vgrad[i]/ve[i]+.5*cfturb(rtheta, h);
    yp[2] = -h1*(vgrad[i]/ve[i]+yp[1]/yt[1])+.0306*pow((h1-3.), -.6169)/yt[1];
}

double h1ofh(double h) {
    if(h > 1.6) {
        return 3.3+1.5501*pow((h-.6778), -3.064); 
    }
    else { 
        return 3.3+.8234*pow((h-1.1), -1.287);
    }
} 

double hofh1(double h1) {
    if(h1 < 5.3) {
        if(h1 < 3.3) {
            return 3.0;
        }
        else {
            return .6778+1.1536*pow((h1-3.3), -.326);
        }
    }
    else {
        return 1.1+.86*pow((h1-3.3), -.777);
    }
}

double cfturb(double rtheta, double h) {
    return .246*pow(10., (-.678*h))*pow(rtheta, -.268);
}