/********************************************************************** 
 * Filename: thinair.ch                                               *
 * Description: This program is coverted from a fortran program       * 
 *    thinair.f used in EME127.  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/18/2007, Debugged and modified by:                          
     Yu-Cheng Chou cycchou@ucdavis.edu                            
  *) 6/4/07,  Converted by 
     Noa Bruhis nbruhis@ucdavis.edu 
**********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <numeric.h>
#include <math.h>
#include <chplot.h>
#define ixx 101
#define TRUE 1
#define EPS     (1.0e-10)
#define PI      (2.0*asin(1.0))
#define DEGRAD  (PI/180.0)

double camberdis(double cxm, double dm, double axii, double teti);
double thickdis(int ithick, double cxm, double em, double axii, double teti);
void moran(int ix, int *itrans, int *isep, double Re, array double x[1:ixx], 
           array double y[1:ixx], array double ve[1:ixx], 
           array double xx[1:ixx], array double vgrad[1:ixx], 
           array double theta[1:ixx], array double cfr[1:ixx]);
void thwats(double alam, double *h, double *el); 
double h1ofh(double h); 
double hofh1(double h1);
double cfturb(double rtheta, double h); 
void derivs(int ix, int i, double Re, array double yp[1:2], 
            array double yt[1:2], array double ve[1:ixx], 
            array double vgrad[1:ixx]);
void runge2(int ix, int i, double Re, double dx, array double yy[1:2], 
            array double ve[1:ixx], array double vgrad[1:ixx]);

int main() {
    double cxm = 1.0, Uref = 1.0, xi, ix, omega;
    double dm, em, ithick, alphd, alpha, dtet;
    int ipr, iter, itrans, isep;
    double teti, axii, fim, eim, axiim, fi;
    int it, itx, idx, i;
    double dgx;
    double sum;
    double Re;
    class CPlot plot[8];
        
    double cpusdi, cposdi;

    array double x[1:ixx], axi[1:ixx];
    array double f[1:ixx], fp[1:ixx];
    array double e[1:ixx], ep[1:ixx];
    array double g[1:ixx], gp[1:ixx], dg[1:ixx];
    array double zu[1:ixx], zo[1:ixx];
    array double ga[1:ixx], gpa[1:ixx];
    array double ue[1:ixx], uf[1:ixx], we[1:ixx], wf[1:ixx];
    array double uu[1:ixx], uo[1:ixx];
    array double wu[1:ixx], wo[1:ixx];
    array double uru[1:ixx], wru[1:ixx];
    array double uro[1:ixx], wro[1:ixx];
    array double cpu[1:ixx], cpo[1:ixx];
    array double cpru[1:ixx], cpro[1:ixx];
    array double ve[1:ixx], xx[1:ixx], vgrad[1:ixx], theta[1:ixx], cfr[1:ixx];

    FILE *fp1, *fp2, *fp3, *fp4, *fp5, *fp6, *fp7, *fp8, *fp9;
    fp1 = fopen("thinair.data", "r");
    fp2 = fopen("thinair.xz", "w");
    fp3 = fopen("thinair.xcp", "w");
    fp4 = fopen("thinair.xu", "w");
    fp5 = fopen("thinair.xw", "w");
    fp6 = fopen("thinair.xg", "w");
    fp7 = fopen("thinair.xgp", "w");
    fp8 = fopen("thinair.xcpsd", "w");

//*****geometric/flow coefficients

    fscanf(fp1, "%lf%*s%*s%*s%*s", &ix);	//%d?
    fscanf(fp1, "%lf%*s%*s%*s", &omega);
    fscanf(fp1, "%lf%*s%*s%*s", &dm);
    fscanf(fp1, "%lf%*s%*s%*s", &em);
    fscanf(fp1, "%lf%*s%*s%*s", &ithick);	//%d?
    fscanf(fp1, "%lf%*s%*s%*s%*s%*s", &alphd);

    alpha = DEGRAD*alphd;
    dtet = PI/(ix-1);

//*****output on listing
    printf("*************************\n");
    printf("dimensionless parameters:\n");
    printf("                           X=C*x\n");
    printf("                           Z=C*z\n");
    printf("                      P-Pinf=0.5*Rho*pow(V,2)*Cp\n");
    printf("                           L=0.5*Rho*pow(V,2)*C*Cl\n");
    printf("                           D=0.5*Rho*pow(V,2)*C*Cd\n");
    printf("                         M,o=0.5*Rho*pow(V,2)*pow(C,2)*Cm,o\n");
    printf("\n");
    printf("******************\n");
    printf("main profile data:\n");
    printf("              mesh system ix=%f\n", ix);
    printf("     relaxation factor omega=%f\n", omega);
    printf("               maximum chord=%f(m)\n", cxm);
    printf("             relative camber=%f\n", dm);
    printf("          relative thickness=%f\n", em);
    printf("      thickness distribution=%f (ellip/semicubic/q-j/naca00xx/selig)\n", ithick);
    printf("***********************\n");
    printf("aerodynamic parameters:\n");
    printf("          reference velocity=%f (m/s)\n", Uref);
    printf("             angle of attack=%f (rd) = %f (deg)\n", alpha, alphd);
    printf(" \n");

//*****initialization
    printf("\n");
    printf("do you want to print the mesh, geometry, and flow? y/n=1/0\n");
    scanf("%d", &ipr);
    if(ipr == 1)
        printf("  i=     x[i]=      zu[i]=       zo[i]=       fp[i]=       ep[i]=\n");
    teti  = 0.5*dtet;
    axii  = 0.5*cxm*(1.0-cos(0.5*dtet));
    fim = camberdis(cxm,dm,axii,teti);
    fi = camberdis(cxm,dm,0,0);
    fim   = 2.0*fi-fim;
    eim = thickdis(ithick,cxm,em,axii,teti);
    axiim = -axii;
    eim   = -eim;

    for(i=1; i<=ix; i++) {
        teti   = (i-1)*dtet;
        xi     = 0.5*cxm*(1.0-cos(teti));
        x[i]   = xi;
        axii   = 0.5*cxm*(1.0-cos(teti+0.5*dtet));
        axi[i] = axii;
        f[i] = camberdis(cxm,dm,axii,teti);
        fp[i]  = (f[i]-fim)/(axii-axiim);

        if(i != ix)
            e[i] = thickdis(ithick,cxm,em,axii,teti+0.5*dtet);
        else
            e[i] = thickdis(ithick,cxm,em,axii,teti-0.5*dtet);

        ep[i] = (e[i] - eim)/(axii - axiim);
        eim   = e[i];
        axii  = 0.5*cxm*(1.-cos(teti));
        e[i] = thickdis(ithick,cxm,em,axii,teti);
        fi = camberdis(cxm,dm,axii,teti);
        zu[i] = fi + e[i] + alpha*(cxm-xi);
        zo[i] = fi - e[i] + alpha*(cxm-xi);
        fim   = f[i];
        f[i]  = fi;
        axiim = axi[i];

        if(i == (ix-1)) {
            fim   = 2.0*fim;
            eim   = 2.0*eim;
            axiim = 2.0*axi[i] - cxm;
        }

        gp[i] = 0;
        g[i]  = 0;

        if(ipr == 1) {
            printf("%3d  %+.4e  %+.4e  %+.4e  %+.4e  %+.4e\n", 
                    i, x[i], zu[i], zo[i], fp[i], ep[i]);
        }
        fprintf(fp2, "%f  %f  %f\n", x[i], zu[i], zo[i]);
    }
    fclose(fp2);

    iter = 0;
    int j;
    double coefi, Cle, dum, sq1, sq2, coefj;
 
//*****iteration loop
    while (TRUE) {
        printf("itx = ?\n");
        scanf("%d", &itx);
        if(itx <= 0)
            break;

        for(it=1; it<=itx; it++) {
            iter += 1;
            idx = 0;
            dgx = 0;

            for(i=2; i<=(ix-1); i++) {
                sum = 0;
                for(j=1; j<=(ix-1); j++) {
                    coefj = 1.0/(x[i] - axi[j]);
                    sum += (g[j+1] - g[j])*coefj;

                    if(j==(i-1))
                         coefi = coefj;
                    if(j==i)
                         coefi -= coefj;
                }
                wf[i] = -sum/(2.0*PI);
                coefi = coefi/(2.0*PI);
                dg[i] = omega*(wf[i] - Uref*(fp[i] - alpha))/coefi;
                g[i] += dg[i];
                if(abs(dgx) < abs(dg[i])) {
                    dgx=dg[i];
                    idx=i;
                }
            }
            g[1]  = 0;
            coefi = pow((x[ix] - x[ix-1]),2)/pow((x[ix] - x[ix-2]),2);
            g[ix] = (coefi*g[ix-2] - g[ix-1])/(coefi - 1.0);
        }
        printf("iter=%d dgx=%f idx=%d\n",iter,dgx,idx);
    }

//*****compute velocity due to camber and thickness
    uf[1] = 0;
    uf[ix] = 0;

    for(i=2; i<=(ix-1); i++) {
        uf[i] = 0.5*(g[i+1] - g[i-1])/(x[i+1] - x[i-1]);
        sum = 0;

        for(j=1; j<=(ix-1); j++) {
            coefj = 1.0/(x[i]-axi[j]);
            sum   += (e[j+1]-e[j])*coefj;
        }
        ue[i] = Uref*sum/PI;
        we[i] = Uref*ep[i];
    }
    we[1] = Uref*ep[1];
    we[ix]= Uref*ep[ix];
//*****exact solution for flat plate at 1.0 rd incidence
    Cle  = 2.0*PI*(1.0+0*2.0*dm/cxm);
    sum  = 0;
    dum  = 0;
    for (i=1; i<=ix; i++) {
        if(i >= 2)
            sum += (uf[i-1]+uf[i])*(x[i]-x[i-1]);
        teti = (i-1)*dtet;
        sq1  = sqrt(axi[i]/cxm);
        sq2  = sqrt(1.0-axi[i]/cxm);
        gpa[i] = 2.0*Uref*(1.0*sq2/sq1+0*8.0*dm*sq1*sq2);

        if(i >= 2)
            dum += 0.5*(gpa[i-1]+gpa[i])*(x[i]-x[i-1]);

        sq1 = sqrt(x[i]/cxm);
        sq2 = sqrt(1.0-x[i]/cxm);
        ga[i] = Uref*cxm*(1.0*(teti+2.0*sq1*sq2)
                +0*2.0*dm*(teti-2.0*sq1*sq2*(1.0-2.0*x[i]/cxm)));   // ?? *0.
    }
    dum  -= 0.5*gpa[ix]*(x[ix] - x[ix-1]);
    coefi = (g[ix] - sum)/dum;
//*****end of iteration-begin post processing
    ue[1] = ue[2]+(ue[3]-ue[2])*(x[1]-x[2])/(x[3]-x[2]);
    ue[ix]= ue[ix-1]+(ue[ix-2]-ue[ix-1])*(x[ix]-x[ix-1])/(x[ix-2]-x[ix-1]);
    wf[1] = Uref*fp[1];
    wf[ix]= Uref*fp[ix];
    for(i=1; i<=(ix-1); i++) {
        gp[i] = (g[i+1]-g[i])/(x[i+1]-x[i]);
        uf[i] += 0.5*coefi*gpa[i];
        uu[i] = ue[i]+uf[i];
        uo[i] = ue[i]-uf[i];
        cpusdi=2.0*uu[i];
        cposdi=2.0*uo[i];
        fprintf(fp8, "%f, %f, %f\n", x[i],cpusdi,cposdi);
        }
    uf[ix] = 0;
    uu[ix] = ue[ix]+uf[ix];
    uo[ix] = ue[ix]-uf[ix];
    gp[ix] = 0;
    cpusdi=2.0*uu[ix];
    cposdi=2.0*uo[ix];
    fprintf(fp8, "%f, %f, %f\n", x[ix],cpusdi,cposdi);
    fclose(fp8);
//*****regularize velocity field due to thickness
    double dx, dz;
    axiim = 0;
    for(i=1; i<=ix; i++) {
        cpu[i] = 2.0*(ue[i]+uf[i])/Uref;
        cpo[i] = 2.0*(ue[i]-uf[i])/Uref;
        dx     = axi[i]-axiim;
        dz     = ep[i]*dx;
        uru[i] = (Uref+ue[i]+uf[i])*dx*dx/(dx*dx+dz*dz);
        wru[i] = uru[i]*dz/dx;
        uru[i] = uru[i]-Uref;
        dz     = -ep[i]*dx;
        uro[i] = (Uref+ue[i]-uf[i])*dx*dx/(dx*dx+dz*dz);
        wro[i] = uro[i]*dz/dx;
        uro[i] = uro[i]-Uref;
        cpru[i] = 2.0*uru[i]/Uref;
        cpro[i] = 2.0*uro[i]/Uref;
        fprintf(fp3, "%f, %f, %f, %f, %f\n", x[i],cpu[i],cpo[i],cpru[i],cpro[i]);
        axiim  = axi[i];

        if(i == ix-1)
            axiim=2.0*axiim-cxm;
    }
    fclose(fp3);

    double Cl, Cmo, Cd, Cdt;
//*****output on listing
    Cl = 2.0*g[ix]/(Uref*cxm);
    Cmo = -0.5*PI*(alpha+2.0*dm);
    Cd = 0;
    coefi = coefi/DEGRAD;
    printf("\n");
    printf("**************inviscid aerodynamic coefficients:\n");
    printf("Cl = %f   Cmo = %f   Cd = %f\n", Cl, Cmo, Cd);
//  printf("gp[i] = %f\n", gp);
//  printf("g[i]=%f\n", g);
//*****output for nxyplot
    for(i=1; i<=ix; i++) {
        fprintf(fp4, "%f, %f, %f, %f, %f\n",x[i],uru[i],uro[i],uu[i],uo[i]);
    }
    fclose(fp4);

//*****output for nxyplot
    for(i=1; i<=ix; i++) {
        wu[i]=Uref*(fp[i]+ep[i]-alpha);
        wo[i]=Uref*(fp[i]-ep[i]-alpha);
        fprintf(fp5, "%f, %f, %f, %f, %f\n", x[i],wru[i],wro[i],wu[i],wo[i]);
        fprintf(fp6, "%f, %f\n",x[i],g[i]);
        fprintf(fp7, "%f, %f\n", x[i],gp[i]);
    }
    fclose(fp5);
    fclose(fp6);
    fclose(fp7);

    //----------------------------------------------------- 
    // Output graphs
    plot[0].data2D(x, zu);
    plot[0].data2D(x, zo);
    plot[0].title("zu(x) and zo(x)");
    plot[0].label(PLOT_AXIS_X, "x");
    plot[0].label(PLOT_AXIS_Y, "zu(x) and zo(x)");
    plot[0].legend("zu(x)", 0);
    plot[0].legend("zo(x)", 1);
    plot[0].plotting();

    plot[1].data2D(x, cpu);
    plot[1].data2D(x, cpo);
    plot[1].title("Cpu(x) and Cpo(x)");
    plot[1].label(PLOT_AXIS_X, "x");
    plot[1].label(PLOT_AXIS_Y, "Cpu(x) and Cpo(x)");
    plot[1].legend("Cpu(x)", 0);
    plot[1].legend("Cpo(x)", 1);
    plot[1].plotting();

    plot[2].data2D(x, uru);
    plot[2].data2D(x, uro);
    plot[2].title("uu(x) and uo(x)");
    plot[2].label(PLOT_AXIS_X, "x");
    plot[2].label(PLOT_AXIS_Y, "uu(x) and uo(x)");
    plot[2].legend("uu(x)", 0);
    plot[2].legend("uo(x)", 1);
    plot[2].plotting();

    plot[3].data2D(x, wru);
    plot[3].data2D(x, wro);
    plot[3].title("wu(x) and wo(x)");
    plot[3].label(PLOT_AXIS_X, "x");
    plot[3].label(PLOT_AXIS_Y, "wu(x) and wo(x)");
    plot[3].legend("wu(x)", 0);
    plot[3].legend("wo(x)", 1);
    plot[3].plotting();

    plot[4].data2D(x, g);
    plot[4].title("g(x)");
    plot[4].label(PLOT_AXIS_X, "x");
    plot[4].label(PLOT_AXIS_Y, "g(x)");
    plot[4].plotting();

    plot[5].data2D(x, gp);
    plot[5].title("gp(x)");
    plot[5].label(PLOT_AXIS_X, "x");
    plot[5].label(PLOT_AXIS_Y, "gp(x)");
    plot[5].plotting();
    //----------------------------------------------------- 
    
    array double x2[1:ixx], cfr2[1:ixx];

//*****boundary layer correction
    printf("\n");
    printf("          ************************\n");
    printf("**********BOUDARY LAYER CORRECTION************\n");
    printf("          ************************\n");

    while(TRUE) {
        printf("\n");
        printf("Reynolds Number Re=?\n");
        scanf("%lf", &Re);
        if(Re < EPS || Re > 1.0/EPS) {
            printf("******input files:\n");
            printf("thinair.data                   :geom. & aero. para.\n");
            printf("******output files:\n");
            printf("thinair.xz                     :x, zu[x] & zo[x]\n");
            printf("thinair.xcp                    :x, Cpu[x] & Cpo[x]\n");
            printf("thinair.xu                     :x, uru[x] & uro[x]
		, uu[x] & uo[x]\n");
            printf("thinair.xw                     :x, wru[x] & wro[x]
		, wu[x] & wo[x]\n");
            printf("thinair.xg                     :x, g[x]\n");
            printf("thinair.xgp                    :x, gp[x]\n");
            printf("thinair.xcfu                   :x, Cfu[x]\n");
            printf("thinair.xcfo                   :x, Cfo[x]\n");
            printf("thinair.xcpsd                  :x, Cpusd[x] & Cposd[x]\n");
            break;
        }
        else {
            printf("\n");
            printf("upper surface\n");
            fp8 = fopen("thinair.xcfu", "w");

            for(i=1; i<=ix; i++) {
                ve[i] = 1.0+uru[i]/Uref;
            }
            moran(ix,&itrans,&isep,Re,x,zu,ve,xx,vgrad,theta,cfr);
            Cd=0;
            for(i=2; i<=ix; i++) {
                Cd=Cd+0.5*(cfr[i-1]+cfr[i])*(x[i]-x[i-1])/cxm;
                if(i > isep)
                    cfr[i]=0;
                fprintf(fp8, "%f, %f\n", x[i], cfr[i]);
                x2[i-1] = x[i];
                cfr2[i-1] = cfr[i];
            }
            fclose(fp8);
    
            //----------------------------------------------------- 
            // Output graphs
            plot[6].data2D(x2, cfr2);
            plot[6].title("Cfu(x)");
            plot[6].label(PLOT_AXIS_X, "x");
            plot[6].label(PLOT_AXIS_Y, "Cfu(x)");
            plot[6].plotting();
            //----------------------------------------------------- 
             
            printf("upper surface itrans=%d isep=%d Cd=%f\n", itrans, isep, Cd);
            Cdt = Cd;
            printf(" \n");
            printf(" \n");
            printf("lower surface\n");
            fp9 = fopen("thinair.xcfo", "w");

            for(i=1; i<=ix; i++)
                ve[i] = 1.0+uro[i]/Uref;

            moran(ix,&itrans,&isep,Re,x,zo,ve,xx,vgrad,theta,cfr);
            Cd=0;

            for(i=2; i<=ix; i++) {
                Cd = Cd+0.5*(cfr[i-1]+cfr[i])*(x[i]-x[i-1])/cxm;
                if(i > isep)
                    cfr[i]=0;
                fprintf(fp9, "%f, %f\n", x[i],cfr[i]);
                x2[i-1] = x[i];
                cfr2[i-1] = cfr[i];
            }
            fclose(fp9);

            //----------------------------------------------------- 
            // Output graphs
            plot[7].data2D(x2, cfr2);
            plot[7].title("Cfo(x)");
            plot[7].label(PLOT_AXIS_X, "x");
            plot[7].label(PLOT_AXIS_Y, "Cfo(x)");
            plot[7].plotting();
            //----------------------------------------------------- 

            printf("lower surface itrans=%d isep=%d Cd=%f\n", itrans, isep, Cd);
            Cdt += Cd;
            printf("total Cd=%f\n", Cdt);
        }
    }
    //scanf("%f", &press);
    return 0;
}


/************** Functions & Subroutines *********************************/

double camberdis(double cxm, double dm, double axii, double teti) {
     double fim;
     fim = 4.0*dm*axii*(1.0-axii/cxm);
     return fim;
}

/************************************************************************/

double thickdis(int ithick, double cxm, double em, double axii, double teti) {
    double fasc = 0.4592793;
    double faqj = 0.3849002;
    double fasl = 0.3088162;
//
//     1=elliptic distribution
//     2=semi-cubic distribution
//     3=quasi-joukowski distribution
//     4=naca00em distribution
//     5=selig distribution
//
    double eim;
    switch(ithick) {
        case 1 :
            eim=0.5*em*cxm*sin(teti);
          break;
        case 2 :
            eim=fasc*em*cxm*sqrt(1.0+cos(teti))*sin(teti);
          break;
        case 3 :
            eim=faqj*em*cxm*(1.0+cos(teti))*sin(teti);
          break;
        case 4 :
            eim=5.0*em*cxm*(0.2969*sqrt(axii/cxm)
                -0.126*axii/cxm-0.3537*pow((axii/cxm),2)
                +0.2843*pow((axii/cxm),3)-0.1015*pow((axii/cxm),4));
            break;
        case 5 :
            eim=fasl*em*cxm*sin(teti)*pow((1.+cos(teti)),1.5);
            break;
        default :
            break;
    }
    return eim;
}

/************************************************************************/

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

    *itrans = ix;
    *isep = ix;
    pi = 3.1415926535;
    eps = 1.0e-12;
    xx[1] = 0;

    for(i=2; i<=ix; 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<=ix; i++) {
        v3 = v1;
        x3 = x1;
        v1 = v2;
        x1 = x2;

        if(i < ix) {
            v2 = ve[i+1];
            x2 = xx[i+1];
        }
        else {
            v2 = ve[ix-1];
            x2 = xx[ix-1];
            v3 = ve[ix-2];
            x3 = xx[ix-2];
        }

        fact     = (x3-x1)/(x2-x1);
        vgrad[i] = ((v2-v1)*fact-(v3-v1)/fact)/(x3-x2);
    }

    printf("     x         y           ve           vdot        theta           h           cf\n");

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

        dv = eps*vgrad[i];
        if(dv < eps*eps)
            dv = eps*eps;

        theta[i] = sqrt(eps*(6.0*pow(v0,5) + 15.0*pow(v0,4)*dv 
                   +20.0*pow(v0,3)*dv*dv+15.0*v0*v0*pow(dv,3)
                   +6.0*v0*pow(dv,4)+pow(dv,5))/(6.0*pow((v0+dv),6)));
        do {
            alam = theta[i]*theta[i]*vgrad[i]*Re;
            if(alam > 0.1)
                alam = 0.1;
            if(alam < -0.0899)
                alam = -0.0899; 
            
            thwats(alam, &h, &el);
 
            if(ve[i] < eps)
                cf = 2.0*el/Re/theta[i]/eps;
            else
                cf = 2.0*el/Re/theta[i]/ve[i];
    
            cfr[i] = cf;
            printf("l %+f, %+f, %+.4e, %+.4e, %+.4e, %+.4e, %+.4e\n", 
                    x[i], y[i], ve[i], vgrad[i], theta[i], h, cf);
            i++;

            if(i > ix)
                return;

            if(ve[i] < 0)
                break;

            dth2ve6  = 0.225*(pow(ve[i],5) + pow(ve[i-1],5))*(xx[i]-xx[i-1])/Re;
            theta[i] = sqrt(((theta[i-1]*theta[i-1])*(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.0+22400.0/Rex)*pow(Rex,0.46);

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

        if(ve[i] < 0) {
            continue;
        }
        else {
            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(ix,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;
        printf("t %f, %f, %+.4e, %+.4e, %+.4e, %+.4e, %+.4e\n", x[i], y[i], ve[i], 
               vgrad[i], theta[i], h, cf);

        if(h > 2.4) {
            *isep = i;
            printf("turbulent separation\n");
            break;
        }

        i++;
    } while(i <= ix);
 
    return; 
}

/************************************************************************/

void thwats(double alam, double *h, double *el) {
    if(alam < 0) {
        *el = 0.22+1.402*alam+0.018*alam/(0.107+alam);
        *h = 2.088+0.0731/(0.14+alam);
    }
    else {
        *el = 0.22+alam*(1.57-1.8*alam);
        *h = 2.61-alam*(3.75-5.24*alam);
    }
}

/************************************************************************/

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

/************************************************************************/

double hofh1(double h1) {
    if(h1 < 5.3) {
        if(h1 < 3.3)
            return 3.0;
        else
            return (pow(0.6778+1.1536*(h1-3.3), (-0.326)));
    }
    else
      return (pow(1.1+0.86*(h1-3.3), (-0.777)));
}

/************************************************************************/

double cfturb(double rtheta, double h) {
    return (pow(0.246*(pow(10.0,(-0.678*h)))*(rtheta), (-0.268)));
}

/************************************************************************/

void derivs(int ix, int i, double Re, array double yp[1:2], array double yt[1:2], 
            array double ve[1:ixx], array double vgrad[1:ixx]) {
    double h1, h, rtheta;
 
    h1   = h1ofh(hofh1(yt[2]));
    h    = hofh1(h1);
    rtheta = Re*ve[i]*yt[1];
    yp[1]  =-(h+2.0)*yt[1]*vgrad[i]/ve[i] + 0.5*cfturb(rtheta,h);
    yp[2]  =pow(-h1*(vgrad[i]/ve[i]+yp[1]/yt[1])+0.0306*(h1-3.0), (-0.6169)/yt[1]);
}

/************************************************************************/

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

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

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

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

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

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