if(%t)
    fd=mopen('constellations_tabbed.txt','r')
    clines=mgetl(fd)
    mclose(fd)
end

dist_eye=640;//[mm] distance from eye to piece of paper
paper_w=11*25.4;//[mm] paper width
paper_h=(8.5+(1/16))*25.4;//[mm] paper height
ang_w=2*atan((paper_w/2)/dist_eye);//[rad] visual angle of the paper
ang_h=2*atan((paper_h/2)/dist_eye);//[rad]

//euler angles defining view field
//for Cassiopeia, -76,-29.5,0 [deg]
//for Ursae Minoris, -40,9,-130 [deg]
//for Ursae Majoris, -83.5,33.8,-58 [deg]
//for Corona Borealis, -34,61,0 [deg]
//for Cygnus, 35,53,0 [deg]
//for Aquila, 22.5,84,14 [deg]
//for Great Square, 79.9,67.5,11.8 [deg]; reduced 84.5,68,3 [deg]
ea=84.5/180*%pi;
eb=68/180*%pi;
ec=3/180*%pi;

//from https://en.wikipedia.org/wiki/Euler_angles
rmat=[cos(ea)*cos(ec)-cos(eb)*sin(ea)*sin(ec),-cos(ea)*sin(ec)-cos(eb)*cos(ec)*sin(ea),sin(ea)*sin(eb);
      cos(ec)*sin(ea)+cos(ea)*cos(eb)*sin(ec),cos(ea)*cos(eb)*cos(ec)-sin(ea)*sin(ec),-cos(ea)*sin(eb);
      sin(eb)*sin(ec),cos(ec)*sin(eb),cos(eb)];
view_dir=(rmat*[0 0 1]')'; //plane normal
up_dir=(rmat*[0 1 0]')'; //these must be orthogonal and unit length
rt_dir=(rmat*[-1 0 0]')';
view_offset=dist_eye; //offset along plane normal

function str=cndisplay(h)
    str=h.user_data;
endfunction

cplot=29;//plot one constellation at a time
cmap=rainbowcolormap(cplot);
cnames=[];
numlines=size(clines)(1);
cind=0;
cpos=[];
drawlater
clf();
ner=0;
for(i=1:numlines)
    st1=clines(i);
    if(isempty(st1))
        continue;
    end
    if(~isdigit(part(st1,1)))
        //this is a constellation name
        if(cind>=cplot)
            break;
        end
        if(%f && (cind>0.5))
            cpos=[mean(cpos,1);cpos];//point for legend entry
            sh=scatter(cpos(:,1),cpos(:,2),36,"markerEdgeColor",cmap(cind,:),"markerFaceColor",cmap(cind,:));
            dh=datatipCreate(sh,1);
            dh.user_data=cnames($);
            datatipSetDisplay(dh,cndisplay);
        end
        if(ner>0.5)
            disp(strcat(["Error: negative projection, ",cnames($)]));
            ner=0;
        end
        cnames($+1)=st1;
        cind=cind+1;
        cpos=[];
    else
        //this is a star
        sparts=strsplit(st1,["	"]);
        ra_h=strtod(sparts(4));
        ra_m=strtod(sparts(5));
        decl=strtod(sparts(6));
        vl=1;
        ra=(ra_h+ra_m/60)/24*(2*%pi);
        dc=decl/360*(2*%pi);
        xs=vl*cos(ra)*cos(dc);
        ys=vl*sin(ra)*cos(dc);
        zs=vl*sin(dc);
        t=view_offset/sum([xs ys zs].*view_dir);
        if(t<0)
            ner=1;
        end
        py=sum((t*[xs ys zs]).*up_dir);
        px=sum((t*[xs ys zs]).*rt_dir);
        cpos($+1,:)=[px py]; //append row
    end
end

cpos=[mean(cpos,1);cpos];//point for legend entry
sh=scatter(cpos(:,1),cpos(:,2),36,"markerEdgeColor",cmap(cind,:),"markerFaceColor",cmap(cind,:));
dh=datatipCreate(sh,1);
dh.user_data=cnames($);
datatipSetDisplay(dh,cndisplay);
if(ner>0.5)
    disp(strcat(["Error: negative projection, ",cnames($)]));
    ner=0;
end

isoview
//gca().data_bounds=[-paper_w/2, -paper_h/2;paper_w/2, paper_h/2]
xrect(-paper_w/2,paper_h/2,paper_w,paper_h);
drawnow
