clear
clc
load -ASCII ah94145.txt
af=ah94145;
k=1:length(af);

top=1:length(af)/2;
bottom=sort(length(af)/2+1:length(af),'descend');
cam_line(:,1)=(af(top,1)+af(bottom,1))/2;
cam_line(:,2)=(af(top,2)+af(bottom,2))/2;
[p o]=max(cam_line(:,2));
max_cam=cam_line(o,:)
% naca=5400;
m=.05;
p=.4;
c=1;

x=linspace(0,1,length(af)/2);
for lw=1:length(af)/2
    for mn=1:length(m)
        for pn=1:length(p)
            if x(lw)<=p
                y_c(lw,mn,pn)=m(mn)/p(pn)^2*(2*p(pn)*x(lw)-x(lw)^2);
            else
                y_c(lw,mn,pn)=m(mn)/(1-p(pn))^2*((1-2*p(pn))+2*p(pn)*x(lw)-x(lw)^2);
            end
        end
    end
end

S = 17.7778;
chord=20/12;
taper=.6;
sweep=deg2rad(.841);
h_span=S./(chord(1,1).*(1+taper(1,:)));
root=chord;
tip=chord*taper;
rib_num=ceil(h_span/(5/12));

% for tx=1:length(taper)
%     x=linspace(0,h_span(tx),rib_num(tx));
%     ribs = (root*(h_span(tx)-x)+tip(tx)*x)./h_span(tx);
%     for rx=1:length(x)
%         area(tx,rx)=polyarea(af(:,1)*ribs(rx),af(:,2)*ribs(rx));
%         [geo iner cpmo]=polygeom(af(:,1)*ribs(rx),af(:,2)*ribs(rx));
%         perim(tx,rx)=geo(4);
%     end
%     wettedarea_b = @(x)interp1(linspace(0,h_span(tx),length(ribs)),perim(tx,:),x);
%     wetted_area(tx) = quad(wettedarea_b,0,h_span(tx));
% end
% disp('     Area          Perimeter')
% disp([area' perim'])
% 
% disp('Wetted Area');disp(wetted_area');
cr_w=2*S./(2*h_span*(1+taper));

chord_x=@(y) cr_w*(taper-1)/(h_span)*y+cr_w;
isodd=@(x) logical(mod(x,2));
rib_chord=chord_x(linspace(0,h_span,rib_num))'*12;
[per_maxthick loc_mt]=max(af(top,2)-af(bottom,2));
for trib=1:rib_num
    rib_para(trib).rib_number=trib;
    rib_para(trib).span_loc=5*(trib-1);
    rib_para(trib).chord=rib_chord(trib);
    rib_para(trib).points=af*rib_chord(trib);
    rib_para(trib).max_thick.size=per_maxthick*rib_chord(trib);
    rib_para(trib).max_thick.chord=af(top(loc_mt),1)*rib_chord(trib);
    rib_para(trib).max_thick.loc_top=af(top(loc_mt),:)*rib_chord(trib);
    rib_para(trib).max_thick.loc_bot=af(bottom(loc_mt),:)*rib_chord(trib);
    rib_para(trib).control_chord=rib_chord(trib)*.80;
    if trib<8
        for num_brk=1
            brk=[.3694 .3694+.5/rib_chord(trib) .80];
            rib_para(trib).rib_break(num_brk).brk=rib_chord(trib)*[brk(num_brk) 2.55/rib_chord(trib);
                brk(num_brk) (2.55-3)/rib_chord(trib);
                brk(num_brk)+.5/rib_chord(trib) 2.55/rib_chord(trib);
                brk(num_brk)+.5/rib_chord(trib) (2.55-3)/rib_chord(trib)];
        end
    elseif trib==8
        for num_brk=1:3
            brk=[.3694 .3694+.5/rib_chord(trib) .80];
            rib_para(trib).rib_break(num_brk).brk=rib_chord(trib)*[brk(num_brk) 2.55/rib_chord(trib);
                brk(num_brk) (2.55-3)/rib_chord(trib);
                brk(num_brk)+.5/rib_chord(trib) 2.55/rib_chord(trib);
                brk(num_brk)+.5/rib_chord(trib) (2.55-3)/rib_chord(trib)];
        end
    elseif trib>8
        for num_brk=1:2
            brk=[.3694+.5/rib_chord(trib) .80];
            rib_para(trib).rib_break(num_brk).brk=rib_chord(trib)*[brk(num_brk) 2.55/rib_chord(trib);
                brk(num_brk) (2.55-3)/rib_chord(trib);
                brk(num_brk)+.5/rib_chord(trib) 2.55/rib_chord(trib);
                brk(num_brk)+.5/rib_chord(trib) (2.55-3)/rib_chord(trib)];
        end
    end
    rb_lc(trib)=rib_para(trib).span_loc;
    szng(trib)=rib_para(trib).max_thick.size;
end
thck=@(x) interp1(1-[1:60]'./60,af(top,2)-af(bottom,2),x);
dspar=@(y) interp1(rb_lc,szng,y);
dh_rad=atan((-dspar(10)+dspar(0))/(10));
dh_ang=rad2deg(dh_rad)

% plot(rib_para(trib).points(:,1),rib_para(trib).points(:,2))
% hold on
% plot(rib_para(trib).points([1;120],1),rib_para(trib).points([1;120],2))
% hold off

% wirte AutoCAD script
wrt=1;
if wrt==1
    start_r=1:rib_num;
    fid = fopen('ribs.scr', 'w');
    quarter_c=[af(:,1)*0 af(:,1)*0 af(:,1)*0];
    line_up=quarter_c;
    dspar_dist=0;
    fprintf(fid, '-units\n2\n2\n1\n1\n\n\n');
    
    if isodd(start_r(1))
    layer_count=0;
        lay_strt=start_r(isodd(start_r));
    else
    layer_count=1;
        lay_strt2=start_r(isodd(start_r));
        lay_strt=[start_r(1) lay_strt2];
            fprintf(fid, '-layer\n');
            fprintf(fid, 'make\n');
            fprintf(fid, 'rib_set%02.0f\n\n',layer_count);
    end
    for trib=start_r
        % Spline the airfoil shape
        rib_para(trib).quarter_c=rib_para(trib).chord*.25;
        rib_para(trib).span_loc=(cr_w*12/4)*(trib-1);
        if trib>1
            dist_c=rib_para(trib-1).chord*.75-rib_para(trib).chord*.75;
            dspar_dist=rib_para(trib-1).max_thick.chord-rib_para(trib).max_thick.chord;
            line_up(:,1)=(trib-1)*dspar_dist+af(:,1)*0;
        end
        if isodd(trib)
            layer_count=layer_count+1;
            line_up(:,2)=0.25+rib_para(1).max_thick.size+af(:,2)*0;
            fprintf(fid, '-layer\n');
            fprintf(fid, 'make\n');
            fprintf(fid, 'rib_set%02.0f\n\n',layer_count);
        else
            line_up(:,2)=af(:,2)*0;
        end
        txt_rib=[line_up(1,:)+[rib_para(trib).max_thick.chord+.5 rib_para(trib).max_thick.size/2,(cr_w*12/4)*(trib-1)];
            line_up(1,:)+[rib_para(trib).control_chord-2 rib_para(trib).max_thick.size*.1,(cr_w*12/4)*(trib-1)];
            line_up(1,:)+[rib_para(trib).control_chord+.25 rib_para(trib).max_thick.size*.1,(cr_w*12/4)*(trib-1)];
            line_up(1,:)+[rib_para(trib).quarter_c-.5 rib_para(trib).max_thick.size/2,(cr_w*12/4)*(trib-1)]];
        fprintf(fid, '-text\n%6.6f,%6.6f,%6.6f\n',txt_rib(1,:));
        fprintf(fid, '.3\n0\nrib%02.0f\n',trib);
        fprintf(fid, '-text\n%6.6f,%6.6f,%6.6f\n',txt_rib(2,:));
        fprintf(fid, '.3\n0\nrib%02.0f\n',trib);
        fprintf(fid, '-text\n%6.6f,%6.6f,%6.6f\n',txt_rib(3,:));
        fprintf(fid, '.15\n0\nrib%02.0f\n',trib);
        fprintf(fid, '-text\n%6.6f,%6.6f,%6.6f\n',txt_rib(4,:));
        fprintf(fid, '.3\n0\nrib%02.0f\n',trib);
        fprintf(fid, 'spline\n');
        fprintf(fid, '%6.6f,%6.6f,%6.6f\n', line_up'+[rib_para(trib).points' ;af(:,1)'*0+(cr_w*12/4)*(trib-1)]);
        fprintf(fid, '\n\n\n');
        r_line=line_up([1 120],:)+[rib_para(trib).points([1 120],1),rib_para(trib).points([1 120],2),af([1 120],1)*0+(cr_w*12/4)*(trib-1)];
        % connect the trailing edge lines with a straight line
        fprintf(fid, 'line\n');
        fprintf(fid, '%6.6f,%6.6f,%6.6f\n', r_line');
        fprintf(fid, '\n');

        % build the rib break lines
        for num_brk=1:length(rib_para(trib).rib_break)
            spar1_line=line_up(1:4,:)+[rib_para(trib).rib_break(num_brk).brk(:,1),rib_para(trib).rib_break(num_brk).brk(:,2), ones(4,1)*rib_para(trib).span_loc];

            fprintf(fid, 'line\n%6.6f,%6.6f,%6.6f\n%6.6f,%6.6f,%6.6f\n\n', spar1_line');
        end
    end
    lay_on=1;
    if layer_count>1
        fprintf(fid, '-layer\non\nrib_set%02.0f\n\n',lay_on);
        fprintf(fid, '-layer\nset\nrib_set%02.0f\n\n',lay_on);
        fprintf(fid, '-layer\noff\nrib_set%02.0f\n\n',2:layer_count);
    end
% Print automatically from AutoCAD    
    prnt=1;
    printer='548 Baldwin manual.pc3';
    if prnt==1
    clear line_up
    
        for trib=1:layer_count
            line_up=[(lay_strt(trib)-1)*dspar_dist 0.25+rib_para(1).max_thick.size];
            if trib>1
                fprintf(fid, '-layer\non\nrib_set%02.0f\n\n',trib);
                fprintf(fid, '-layer\nset\nrib_set%02.0f\n\n',trib);
                fprintf(fid, '-layer\noff\nrib_set%02.0f\n\n',lay_on);
                lay_on=trib;
            end
            %printing the front part of the ribs from AutoCAD
            fprintf(fid, '-plot\ny\n\n%s\nLetter\nInches\nLandscape\nNo\nWindow\n',printer);
            fprintf(fid, '%6.6f,%6.6f\n',[-.1 min(rib_para(1).points(:,2)-.1);
                line_up(1)+rib_para(lay_strt(trib)).rib_break(1).brk(3,1) line_up(2)+rib_para(lay_strt(trib)).rib_break(1).brk(3,2)]');
            fprintf(fid, '1=1\n\n\n\n\n\n\n\n\n');

            %printing the rear part of the ribs from AutoCAD
            fprintf(fid, '-plot\ny\n\n%s\nLetter\nInches\nLandscape\nNo\nWindow\n',printer);
            fprintf(fid, '%6.6f,%6.6f\n',[line_up(1)+rib_para(lay_strt(trib)).rib_break(1).brk(3,1) min(rib_para(1).points(:,2)-.1);
                rib_para(lay_strt(trib)).chord+line_up(1)+.1 line_up(2)+max(rib_para(1).points(:,2))]');
            fprintf(fid, '1=1\n\n\n\n\n\n\n\n\n');
            
            %printing the very rear part of the ribs from AutoCAD
            if trib>0&trib<4
            fprintf(fid, '-plot\ny\n\n%s\nLetter\nInches\nLandscape\nNo\nWindow\n',printer);
            fprintf(fid, '%6.6f,%6.6f\n',[line_up(1)+rib_para(lay_strt(trib)).rib_break(1).brk(3,1) min(rib_para(1).points(:,2)-.1);
                rib_para(lay_strt(trib)).chord+line_up(1)+.1 line_up(2)+max(rib_para(1).points(:,2))]'+[4 4;0 0]);
            fprintf(fid, '1=1\n\n\n\n\n\n\n\n\n');
            end
        end
    end

    fclose(fid);
end