program specific_speed
    implicit none
    integer::mtype,iostat_val,i,n_pts
    double precision::Q,H,N,D,omega,Ns,Ds,Ns_US,eta_cord,Ds_cord,D_cord
    double precision::g,rho_air,H_m,N_sw,Ns_sw,pi_v
    character(len=40)::mname,regime
    pi_v=3.14159265358979d0;g=9.81d0;rho_air=1.2d0
    read(*,*,iostat=iostat_val) mtype;read(*,*,iostat=iostat_val) Q
    read(*,*,iostat=iostat_val) H;read(*,*,iostat=iostat_val) N
    read(*,*,iostat=iostat_val) D
    if(iostat_val/=0)then;write(*,*)'ERROR: Bad input.';stop;end if
    if(Q<=0)Q=0.1d0;if(H<=0)H=10.0d0;if(N<=0)N=1450.0d0
    select case(mtype)
    case(1);mname='Centrifugal Pump';H_m=H
    case(2);mname='Axial Pump';H_m=H
    case(3);mname='Centrifugal Fan';H_m=H/(rho_air*g)
    case(4);mname='Axial Fan';H_m=H/(rho_air*g)
    case(5);mname='Radial Turbine';H_m=H
    case default;mname='Axial Turbine';H_m=H;mtype=6
    end select
    omega=N*2.0d0*pi_v/60.0d0
    Ns=omega*sqrt(Q)/(g*H_m)**0.75d0
    Ns_US=Ns*2733.0d0
    if(D>0)then
      Ds=D*(g*H_m)**0.25d0/sqrt(Q)
    else
      Ds_cord=exp(1.25d0-0.40d0*log(max(Ns,0.01d0)))
      Ds=Ds_cord
      D_cord=Ds*sqrt(Q)/(g*H_m)**0.25d0
      D=D_cord
    end if
    Ds_cord=exp(1.25d0-0.40d0*log(max(Ns,0.01d0)))
    eta_cord=0.90d0-0.08d0*(log(Ns)-log(1.0d0))**2
    if(eta_cord<0.40d0)eta_cord=0.40d0;if(eta_cord>0.95d0)eta_cord=0.95d0
    if(Ns<0.2d0)then;regime='Positive displacement'
    else if(Ns<1.0d0)then;regime='Radial / Centrifugal'
    else if(Ns<4.0d0)then;regime='Mixed flow'
    else;regime='Axial flow';end if
    write(*,'(A)')'============================================================'
    write(*,'(A)')'   SPECIFIC SPEED / SPECIFIC DIAMETER (Cordier)'
    write(*,'(A)')'============================================================'
    write(*,*)
    write(*,'(A)')'--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')       '  Machine Type              = ',trim(mname)
    write(*,'(A,F10.4,A)') '  Volume Flow Q             = ',Q,' m3/s'
    write(*,'(A,F10.2,A)') '  Head H                    = ',H,' m (or Pa for fans)'
    write(*,'(A,F10.2,A)') '  Head (equiv) H_m          = ',H_m,' m'
    write(*,'(A,F10.1,A)') '  Speed N                   = ',N,' rpm'
    write(*,'(A,F10.4,A)') '  Diameter D                = ',D,' m'
    write(*,*)
    write(*,'(A)')'--- SPECIFIC SPEED ------------------------------------------'
    write(*,'(A,F12.4)')   '  Ns (dimensionless)        = ',Ns
    write(*,'(A,F12.1)')   '  Ns (US customary)         = ',Ns_US
    write(*,'(A,F12.4)')   '  Ds (specific diameter)    = ',Ds
    write(*,'(A,F12.4)')   '  omega [rad/s]             = ',omega
    write(*,*)
    write(*,'(A)')'--- CORDIER DIAGRAM -----------------------------------------'
    write(*,'(A,F12.4)')   '  Ds_Cordier (predicted)    = ',Ds_cord
    write(*,'(A,F12.4)')   '  eta_Cordier (predicted)   = ',eta_cord
    write(*,'(A,A)')       '  Flow Regime               = ',trim(regime)
    write(*,*)
    n_pts=25
    write(*,'(A)')'--- SWEEP: Ns VS SPEED N ------------------------------------'
    write(*,'(A)')'  N[rpm]      omega       Ns          Ns_US       Regime'
    write(*,'(A)')'  -----------------------------------------------------------'
    do i=1,n_pts
      N_sw=500.0d0+dble(i-1)*(5000.0d0-500.0d0)/dble(n_pts-1)
      omega=N_sw*2.0d0*pi_v/60.0d0
      Ns_sw=omega*sqrt(Q)/(g*H_m)**0.75d0
      if(Ns_sw<0.2d0)then;regime='PD'
      else if(Ns_sw<1.0d0)then;regime='Radial'
      else if(Ns_sw<4.0d0)then;regime='Mixed'
      else;regime='Axial';end if
      write(*,'(F8.0,4X,F8.2,4X,F10.4,4X,F10.1,4X,A)') N_sw,omega,Ns_sw,Ns_sw*2733.0d0,trim(regime)
    end do
    write(*,*)
    write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)')'  Ns = omega*sqrt(Q) / (g*H)^0.75'
    write(*,'(A)')'  Ds = D*(g*H)^0.25 / sqrt(Q)'
    write(*,'(A)')'  Cordier: Ds = exp(1.25 - 0.40*ln(Ns))'
end program specific_speed