/********************************************************************* 
 *
 *  Fan mesh definition
 * 
 *  Set the blade length, number of blade stations, and the airfoil
 *  parameters at each station, and the number of blades, and Gmsh
 *  will create the O-grid fan block. Airfoil coordinates are defined
 *  in airfoil.geo.
 * 
 *********************************************************************/

Include "airfoil.geo";

// Fuselage Station of the Fan
FanStation = 42;

// Blade length in inches
BladeLength = 11;

// Point along chord where the axial blade reference curve passes
// through the plane of the airfoil. Given as a percentage of chord
// length.
AirfoilCenter = 0.3;

// Hub radius in inches
HubForeRadius = 4;
HubAftRadius = 3;
HubThickness = 4;

// Number of blades
BladeCount = 7;

// Number of blade stations, minimum of two
StationCount = 5;

ChordMesh = 17;
AxialMesh = 9;


//
Chord[1] = 3.5;
Airfoil[1] = MH112;
Angle[1] = 78;

//Chord[2] = 3.375;
//Airfoil[2] = MH112;
//Angle[2] = 74;

Chord[2] = 3.25;
Airfoil[2] = MH113;
Angle[2] = 68;

//Chord[4] = 3.125;
//Airfoil[4] = MH114;
//Angle[4] = 64;

Chord[3] = 3.0;
Airfoil[3] = MH115;
Angle[3] = 59;

//Chord[6] = 2.875;
//Airfoil[6] = MH116;
//Angle[6] = 55;

Chord[4] = 2.75;
Airfoil[4] = MH117;
Angle[4] = 51;

//Chord[8] = 2.625;
//Airfoil[8] = MH120;
//Angle[8] = 47;

Chord[5] = 2.5;
Airfoil[5] = MH121;
Angle[5] = 44;

HubAftCenter = newp;
Point (HubAftCenter) = {FanStation + HubThickness / 2, 0, 0, 1};

HubForeCenter = newp;
Point (HubForeCenter) = {FanStation - HubThickness / 2, 0, 0, 1};

For blade In {0:BladeCount-1}

  // Global rotation angle of the blade local z axis
  omega = 2 * Pi * blade/BladeCount;

  HubAftArcStart[blade] = newp;
  Point (HubAftArcStart[blade]) = {FanStation + HubThickness / 2,
                                   Sin(omega - Pi / BladeCount) * HubAftRadius,
                                   Cos(omega - Pi / BladeCount) * HubAftRadius, 1};

  HubAftArcEnd[blade] = newp;
  Point (HubAftArcEnd[blade]) = {FanStation + HubThickness / 2,
                                   Sin(omega + Pi / BladeCount) * HubAftRadius,
                                   Cos(omega + Pi / BladeCount) * HubAftRadius, 1};

  HubForeArcStart[blade] = newp;
  Point (HubForeArcStart[blade]) = {FanStation - HubThickness / 2,
                                   Sin(omega - Pi / BladeCount) * HubForeRadius,
                                   Cos(omega - Pi / BladeCount) * HubForeRadius, 1};

  HubForeArcEnd[blade] = newp;
  Point (HubForeArcEnd[blade]) = {FanStation - HubThickness / 2,
                                   Sin(omega + Pi / BladeCount) * HubForeRadius,
                                   Cos(omega + Pi / BladeCount) * HubForeRadius, 1};

  HubAftArc[blade] = newl;
  Circle (HubAftArc[blade]) = {HubAftArcStart[blade], HubAftCenter, HubAftArcEnd[blade]};

  HubForeArc[blade] = newl;
  Circle (HubForeArc[blade]) = {HubForeArcStart[blade], HubForeCenter, HubForeArcEnd[blade]};

  HubSplit1[blade] = newl;
  Line (HubSplit1[blade]) = {HubAftArcStart[blade], HubForeArcStart[blade]};

  HubSplit2[blade] = newl;
  Line (HubSplit2[blade]) = {HubAftArcEnd[blade], HubForeArcEnd[blade]};

  HubPatch[blade] = news;
  Line Loop (HubPatch[blade]) = {HubAftArc[blade],
                      HubSplit2[blade],
                      -HubForeArc[blade],
                      -HubSplit1[blade]};

  Ruled Surface (HubPatch[blade] + 1) = {HubPatch[blade]};

  For station In {1:StationCount}

    shape = Airfoil[station];
    index = shape * IndexFactor;

    scale = Chord[station];
    theta = ((Angle[station]+270) / 360) * (2 * Pi);
    radius = (HubForeRadius + HubAftRadius) / 2 + 
             (station - 1) * BladeLength / (StationCount - 1);

    AirfoilTopStartPoint[station] = newp;

    For p In {0:CoordCnt[shape] - 1}

      af_lx = scale * (Coordinate[index+(4*p)] - AirfoilCenter);
      af_ly = scale * Coordinate[index+(4*p)+1];

      rot_lx = af_lx * Cos(theta) + af_ly * Sin(theta);
      rot_ly = af_ly * Cos(theta) - af_lx * Sin(theta);


      // At the root, conform the airfoil to the hub
      If (station == 1)
        radius = HubForeRadius - (HubForeRadius - HubAftRadius) / HubThickness * 
                 (rot_lx - scale * (0.5 - AirfoilCenter) + (HubThickness / 2));
      EndIf

      delta_omega = Atan(rot_ly/radius);

      gx = rot_lx + FanStation - scale * (0.5 - AirfoilCenter);
      gy = Sin(omega + delta_omega) * radius;
      gz = Cos(omega + delta_omega) * radius;

      // Find the leading edge
      If ((Coordinate[index+(4*p)] < Coordinate[index+(4*(p-1))]) &&
          (Coordinate[index+(4*p)] < Coordinate[index+(4*(p+1))]))
        AirfoilTopEndPoint[station] = newp;
        AirfoilBottomStartPoint[station] = newp;
      EndIf

      Point(newp) = {gx,
                     gy,
                     gz,
                     Coordinate[index+4*p+3]};

    EndFor

    AirfoilBottomEndPoint[station] = newp - 1;

    AirfoilTopSpline[station] = newl;
    Spline(AirfoilTopSpline[station]) = {AirfoilTopStartPoint[station]:AirfoilTopEndPoint[station]};
    Transfinite Line {AirfoilTopSpline[station]} = ChordMesh Using Progression .85;

    AirfoilBottomSpline[station] = newl;
    Spline(AirfoilBottomSpline[station]) = {AirfoilBottomStartPoint[station]:AirfoilBottomEndPoint[station]};
    Transfinite Line {AirfoilBottomSpline[station]} = ChordMesh Using Progression 1.2;

  EndFor
  
  For station In {2:StationCount}

    TopTrailingEdge[station] = newl;
    Line (TopTrailingEdge[station]) = {AirfoilTopStartPoint[station-1],
                                       AirfoilTopStartPoint[station]};
    Transfinite Line {TopTrailingEdge[station]} = AxialMesh;

    BottomTrailingEdge[station] = newl;
    Line (BottomTrailingEdge[station]) = {AirfoilBottomEndPoint[station-1],
                                          AirfoilBottomEndPoint[station]};
    Transfinite Line {BottomTrailingEdge[station]} = AxialMesh;

    TopLeadingEdge[station] = newl;
    Line (TopLeadingEdge[station]) = {AirfoilTopEndPoint[station-1],
                                      AirfoilTopEndPoint[station]};
    Transfinite Line {TopLeadingEdge[station]} = AxialMesh;

    BottomLeadingEdge[station] = newl;
    Line (BottomLeadingEdge[station]) = {AirfoilBottomStartPoint[station-1],
                                         AirfoilBottomStartPoint[station]};
    Transfinite Line {BottomLeadingEdge[station]} = AxialMesh;

    BladePatch[station] = newreg;
    Line Loop (BladePatch[station]) = {AirfoilTopSpline[station-1],
                                       TopLeadingEdge[station],
                                       -AirfoilTopSpline[station],
                                       -TopTrailingEdge[station]};

    Ruled Surface (BladePatch[station] + 1) = {BladePatch[station]};
    Transfinite Surface {BladePatch[station] + 1} = {AirfoilTopStartPoint[station-1],
                                                     AirfoilTopEndPoint[station-1],
                                                     AirfoilTopEndPoint[station],
                                                     AirfoilTopStartPoint[station]};

    Recombine Surface {BladePatch[station] + 1};

    Line Loop (BladePatch[station] + 2) = {AirfoilBottomSpline[station-1],
                                           BottomTrailingEdge[station],
                                          -AirfoilBottomSpline[station],
                                          -BottomLeadingEdge[station]};

    Ruled Surface (BladePatch[station] + 3) = {BladePatch[station] + 2};
    Transfinite Surface {BladePatch[station] + 3} = {AirfoilBottomStartPoint[station-1],
                                                     AirfoilBottomEndPoint[station-1],
                                                     AirfoilBottomEndPoint[station],
                                                     AirfoilBottomStartPoint[station]};

    Recombine Surface {BladePatch[station] + 3};

  EndFor  

EndFor