using CCBlade
using Statistics: mean

# Comparison using geometry and data from)
# https://rotorcraft.arc.nasa.gov/Publications/files/RamasamyGB_ERF10_836.pdf


# add an "overloaded" function to handle the tip loss used just for hover)
function solvehover(rotor, section, op)

	 if section.r/rotor.Rtip > 0.94
        phi = 0.0
        theta = section.theta + op.pitch
        alpha = theta - phi
        _, cd = afeval(section.af, alpha, 1.0, 1.0)
        cl = 0.0

        # resolve into normal and tangential forces
        cn = 0.0
        ct = cd
        F = 1.0
        G = 1.0

        u = 0.0
        v = 0.0
        a = 0.0
        ap = 0.0
        W = op.Vy
        Np = 0.0
        Tp = ct*0.5*op.rho*W^2*section.chord

        return Outputs(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G)
    else
    	return solve(rotor, section, op)
    end
end


chord = 0.060
theta = 0.0
Rtip = 0.656
Rhub = 0.19*Rtip

# turbine = false
B = 3

rotor = Rotor(Rhub, Rtip, B, turbine=false)

r = range(Rhub, Rtip, length=30)

af = AlphaAF("naca0012v2.txt", radians=false)
function af2(alpha, Re, M)
	cl, cd = afeval(af, alpha, Re, M)
	return cl, cd+0.014
end
sections = Section.(r, chord, theta, af2)


rho = 1.225
Omega = 800*pi/30
Vinf = 0.0

Vinf2 = 1e-4

nP = 40
pitch = range(1e-4, 20*pi/180, length=nP)

CT = zeros(nP)
CQ = zeros(nP)
FM = zeros(nP)
CT2 = zeros(nP)
CQ2 = zeros(nP)
FM2 = zeros(nP)


for i = 1:nP
	local op = simple_op.(Vinf, Omega, r, rho, pitch=pitch[i])
	local op2 = simple_op.(Vinf2, Omega, r, rho, pitch=pitch[i])
	local outputs = solvehover.(Ref(rotor), sections, op)
	local outputs2 = solvehover.(Ref(rotor), sections, op2)
	local T, Q = thrusttorque(rotor, sections, outputs)
	T2, Q2 = thrusttorque(rotor, sections, outputs2)

	FM[i], CT[i], CQ[i] = nondim(T, Q, Vinf, Omega, rho, rotor, "helicopter")
	FM2[i], CT2[i], CQ2[i] = nondim(T2, Q2, Vinf, Omega, rho, rotor, "helicopter")
	
end

sigma = B * chord / (pi * Rtip)


using PyPlot
close("all")

rc("axes.spines", right=false, top=false)
rc("font", size=14.0)
rc("legend", frameon=false)
rc("text", usetex="True")


figure()
plot(CT/sigma, CQ/sigma, color="#348ABD")
# plot(CT2/sigma, CQ2/sigma, "k--")
xlim([0, 0.15])
ylim([0, 0.025])
xlabel(L"C_T/\sigma")
ylabel(L"C_Q/\sigma", rotation=0)

data = [
0.0006486486486486538 0.001826767151767151
0.00021621621621621678 0.0019711018711018706
0.0009729729729729721 0.0024036122661122636
0.0021621621621621644 0.002283134095634097
0.010594594594594595 0.0027378378378378337
0.011351351351351346 0.0026174636174636154
0.012540540540540546 0.0024008316008316023
0.011675675675675679 0.003098154885654885
0.012000000000000004 0.0034346153846153867
0.017081081081081077 0.0035535862785862765
0.01935135135135136 0.0036491943866943863
0.017621621621621626 0.0029765332640332605
0.017081081081081077 0.0028324324324324295
0.019135135135135137 0.0027357848232848203
0.05016216216216216 0.004891787941787941
0.05059459459459458 0.005228222453222447
0.054054054054054064 0.005684121621621619
0.05502702702702702 0.006501195426195424
0.05427027027027025 0.006597531185031182
0.06940540540540538 0.007194854469854465
0.07189189189189188 0.007410602910602908
0.07243243243243241 0.007771049896049893
0.07178378378378378 0.00810774428274428
0.07675675675675675 0.008322895010395012
0.0829189189189189 0.008657952182952181
0.08572972972972973 0.009017853430353426
0.08637837837837836 0.0090176975051975
0.088 0.00923365384615384
0.08637837837837836 0.009570582120582115
0.0937297297297297 0.010602468814968808
0.09448648648648647 0.010265748440748436
0.09848648648648647 0.010841709979209978
0.0963243243243243 0.011178768191268186
0.09448648648648647 0.01264555613305613
0.11589189189189186 0.014948102910602907
0.11794594594594593 0.015139916839916836
0.11816216216216216 0.015284095634095628
0.11805405405405403 0.015596621621621619
0.11502702702702702 0.015621387733887729
0.11956756756756756 0.014586642411642408
0.12054054054054059 0.014850831600831596
0.1211891891891892 0.014850675675675672
]
plot(data[:, 1], data[:, 2], "o", color="#A60628")
text(0.133, 0.014, "BEM", color="#348ABD")
text(0.08, 0.017, "experimental", color="#A60628")
gca().yaxis.set_label_coords(-0.09,1.05)
# savefig("hoverCT.pdf")


figure()
plot(CT/sigma, FM, color="#348ABD")
xlim([0, 0.15])
ylim([0, 0.7])
xlabel(L"C_T/\sigma")
ylabel("figure of merit", rotation=0)

data = [
0.0007194244604316561 0.0031982942430702765
0.011305241521068862 0.09786780383795313
0.017163412127440908 0.16439232409381654
0.05056526207605345 0.45223880597014926
0.07194244604316546 0.5462686567164178
0.08581706063720454 0.5827292110874199
0.09362795477903393 0.5641791044776119
0.12127440904419323 0.5936034115138591
0.11582733812949644 0.5533049040511726
]
plot(data[:, 1], data[:, 2], "o", color="#A60628")
text(0.12, 0.67, "BEM", color="#348ABD")
text(0.12, 0.5, "experimental", color="#A60628")
gca().yaxis.set_label_coords(0,1.05)
# savefig("hoverFM.pdf")



# figure()
# plot(CT/sigma, CQ/sigma, color="#348ABD")
# plot(CT2/sigma, CQ2/sigma, "k--")
# xlim([0, 0.15])
# ylim([0, 0.025])
# xlabel(L"C_T/\sigma")
# ylabel(L"C_Q/\sigma", rotation=0)
# gca().yaxis.set_label_coords(-0.09,1.05)
# legend(["modified", "original"])
# # savefig("hoverCTcompare.pdf")

# figure()
# plot(CT/sigma, FM, color="#348ABD")
# plot(CT2/sigma, FM2, "k--")
# xlim([0, 0.15])
# ylim([0, 0.7])
# xlabel(L"C_T/\sigma")
# ylabel("figure of merit", rotation=0)
# gca().yaxis.set_label_coords(0,1.05)
# legend(["modified", "original"])
# # savefig("hoverFMcompare.pdf")

# figure()
# plot(pitch*180/pi, CT/sigma, color="#348ABD")
# plot(pitch*180/pi, CT2/sigma, "k--")
# xlabel("pitch (deg)")
# ylabel(L"C_T/\sigma", rotation=0)
# gca().yaxis.set_label_coords(-0.09,1.05)

figure()
plot(pitch*180/pi, CQ/sigma, color="#348ABD")
plot(pitch*180/pi, CQ2/sigma, "k--")
xlabel("pitch (deg)")
ylabel(L"C_Q/\sigma", rotation=0)
gca().yaxis.set_label_coords(-0.09,1.05)
text(13.7, 0.025, "reformulated", color="#348ABD")
text(18.5, 0.02, "original")
# savefig("CQcompare.pdf")



# ---------- distributed loads --------


# adjust pitch to match thrust
# CT = 0.0013
op = simple_op.(Vinf, Omega, r, rho, pitch=3.5*pi/180)
outputs = solvehover.(Ref(rotor), sections, op)

T, Q = thrusttorque(rotor, sections, outputs)
A = pi*Rtip^2
CT = T / (rho * A * (Omega * Rtip)^2)
println(CT)

figure()
plot(r/Rtip, outputs.u/(Omega*Rtip), color="#348ABD")
xlabel(L"r/R")
ylabel(L"u/(\Omega R)", rotation=0)

data = [
0.20619554695062925 -0.011297248925009579
0.21055179090029041 -0.009261520970755763
0.21926427879961274 -0.00896913482969
0.22434656340755083 -0.013037917332672955
0.24394966118102612 -0.011289932236205238
0.2323330106485963 -0.014780555618091357
0.23814133591481124 -0.013325941601566901
0.25338818973862537 -0.013759033296562279
0.2642787996127783 -0.015065062248137043
0.27516940948693125 -0.016225742362502543
0.28751210067763794 -0.015641955019248532
0.297676669893514 -0.018256264211260947
0.3136495643756051 -0.0188345640379117
0.3238141335914811 -0.02057678020666831
0.3354307841239109 -0.019121040545712426
0.3463213939980639 -0.020572418334496495
0.36592449177153913 -0.02245815416826133
0.37536302032913843 -0.02129353429838582
0.38625363020329134 -0.022744912087169916
0.39932236205227484 -0.023759821247664292
0.41384317521781216 -0.022884914113329907
0.4283639883833494 -0.023754193025507117
0.4428848015488867 -0.02564091379814945
0.46030977734753137 -0.02578288570206444
0.47773475314617614 -0.027087648303653858
0.49661181026137463 -0.028246780656926096
0.5169409486931267 -0.028388189738625377
0.537270087124879 -0.028965645331952575
0.5546950629235237 -0.030125059096332674
0.573572120038722 -0.030993493775186315
0.5793804453049369 -0.029684228595871137
0.5939012584704744 -0.030553507508048375
0.6084220716360116 -0.029823949210923253
0.62729912875121 -0.02996563970373041
0.6454501452081316 -0.028363284855579846
0.6606969990319458 -0.029377771899412444
0.6803000968054211 -0.029664670523874948
0.70135527589545 -0.02980593890002027
0.7187802516940947 -0.030093259641144593
0.739109390125847 -0.03023466872284386
0.7550822846079381 -0.03153971273554111
0.7681510164569214 -0.03269997073324481
0.776863504356244 -0.03168084040613252
0.7899322362052275 -0.031969005380580395
0.7971926427879961 -0.03298504018550623
0.813165537270087 -0.032400549314482566
0.8349467570183928 -0.03268702582228328
0.8480154888673765 -0.03370193498277767
0.8654404646660212 -0.03384390688669267
0.8741529525653435 -0.033260823071208284
0.8785091965150047 -0.03151579279137311
0.8828654404646661 -0.029334715999909985
0.8850435624394966 -0.02729941016231796
0.887947725072604 -0.0251186147819627
0.8952081316553726 -0.021919533307818748
0.9010164569215875 -0.017993989058736143
0.9060987415295256 -0.014213934352416788
0.9126331074540173 -0.010433598234989544
0.9184414327202324 -0.005781309799860454
0.9242497579864473 -0.000983672527522031
0.9300580832526621 0.0032325693959791638
0.9373184898354305 0.007449092730588226
0.9431268151016456 0.011229288142461526
0.9496611810261375 0.01457357774826084
0.952565343659245 0.01588228010536029
0.9612778315585673 0.01588396857200744
0.9699903194578896 0.016321703550282514
0.9787028073572119 0.015015252482045938
0.9903194578896417 0.012982620049978577
1.0019361084220715 0.01327556901326006
0.9837850919651501 0.013998795560458334
]
plot(data[:, 1], -data[:, 2], "o", color="#A60628")
gca().invert_yaxis()
text(0.5, 0.015, "BEM", color="#348ABD")
text(0.5, 0.042, "experimental", color="#A60628")
xlim([0, 1])
ylim([0.1, -0.05])
gca().yaxis.set_label_coords(-0.08,1.05)
# savefig("hoveru1.pdf")


# CT = 0.0076
op = simple_op.(Vinf, Omega, r, rho, pitch=11.8*pi/180)
outputs = solvehover.(Ref(rotor), sections, op)

T, Q = thrusttorque(rotor, sections, outputs)
A = pi*Rtip^2
CT = T / (rho * A * (Omega * Rtip)^2)
println(CT)

figure()
plot(r/Rtip, outputs.u/(Omega*Rtip), color="#348ABD")
xlabel(L"r/R")
ylabel(L"u/(\Omega R)", rotation=0)

data = [
0.19326598413451757 -0.008918299545731909
0.20112550003390056 -0.011346396365855313
0.20542850362736453 -0.01491436707573393
0.20972642213031384 -0.017769069089429795
0.22470811580446126 -0.01920130178317174
0.2147321852328971 -0.01991077361177028
0.22757610685470198 -0.021484846430266435
0.23614143331751297 -0.022914638280561395
0.23758356498745675 -0.025197640518001233
0.2432809004000271 -0.02434388772120144
0.2525733948064275 -0.027771103125635632
0.25826056003796866 -0.025490812936470264
0.26537968675842427 -0.024066987592379144
0.27110956675028813 -0.02777815445114923
0.28323547359142986 -0.02863868736863516
0.29180385110854967 -0.030496440436639766
0.306067530002034 -0.031215133229371472
0.3146318394467421 -0.032502271340429864
0.3274737270323411 -0.03379103668045291
0.34173130381720784 -0.03365380703776527
0.35244457251339073 -0.03636829615567157
0.3659963387348294 -0.037229371482812385
0.37598955861414324 -0.03894501322123532
0.39026340768865686 -0.04109024340633262
0.4045230185097294 -0.0412383212421181
0.41735270187809337 -0.04081524171130245
0.43020170859041273 -0.043102583225981414
0.44304563021221766 -0.04467665604447758
0.4630178317174045 -0.04611078717201167
0.48156213980608853 -0.04725906841141772
0.4986846565868872 -0.048977422198115114
0.5257800528849412 -0.04955834293850429
0.5514516238389042 -0.05042402874771171
0.577113024611838 -0.04986317716455352
0.5956512305919045 -0.050155535968540235
0.6198918570750558 -0.05030741067190994
0.6441436707573395 -0.05202847650688182
0.6826459420977692 -0.05261373652451014
0.70404298596515 -0.05390575632246254
0.7211563495830225 -0.05434022645603091
0.7439751847582886 -0.055062173706691975
0.7610936334666756 -0.056209912536443155
0.7767852057766627 -0.05721445521730287
0.7981822496440435 -0.05850647501525526
0.8081795375957691 -0.06079273171062444
0.816750966167198 -0.06307844599633873
0.8267574072818495 -0.06664858634483693
0.841019052139128 -0.06708197165909553
0.8524096548918569 -0.06480385110854973
0.8609424367753749 -0.06166872330327479
0.8637738151739102 -0.05881673333785341
0.8680188487355073 -0.054253440911248216
0.8708329378262931 -0.04897633737880533
0.8729239270459013 -0.042272425249169424
0.879988134788799 -0.033145297986304154
0.8849298257508984 -0.02629981693674148
0.8912983931113971 -0.01959753203607023
0.8969519967455417 -0.012609668452098446
0.9033215811241437 -0.006050037290663765
0.9089792528307004 0.0003672113363617918
0.9146409926096684 0.006213845006441114
0.9210166790968879 0.01191755373245644
0.926677401857753 0.01790684114177233
0.9330551223811783 0.023325242389314535
0.9394348769408092 0.028458336158383626
0.9443948742287611 0.032736049901688255
0.9486602481524171 0.03444626754356228
0.9543535154925755 0.035870635297308294
0.9629005356295339 0.037008610753271415
0.9728856193640245 0.03643419892874094
0.9786032951386534 0.03443487694080956
0.9885944809817613 0.03300454268085973
0.9957339480642755 0.03157529324021968
0.998601939114516 0.029291748593124965
]
plot(data[:, 1], -data[:, 2], "o", color="#A60628")
gca().invert_yaxis()
text(0.5, 0.068, "BEM", color="#348ABD")
text(0.5, 0.04, "experimental", color="#A60628")
xlim([0, 1])
ylim([0.1, -0.05])
gca().yaxis.set_label_coords(-0.08,1.05)
# savefig("hoveru2.pdf")