とりあえず、後半部分をこうしてみた。
if((xx1.realp() && para1.realp()) || (!xx1.realp() && !para1.realp())){
if(n1+0.01<n2){
POutput->xx = xx1;
para1 = xx1;
POutput->yy = yy1;
para2 = yy1;
}
else {
POutput->xx = xx2;
para1 = xx2;
POutput->yy = yy2;
para2 = yy2;
}
POutput->zz = CComplex(1.0,0.0);
}
else {
if(norm(xx1-xx2)>norm(yy1-yy2)){
px1=(xx1-para1)/(xx2-para1);
}
else {
px1=(yy1-para2)/(yy2-para2);
}
if(px1.imag()<0){
POutput->xx = xx1;
para1 = xx1;
POutput->yy = yy1;
para2 = yy1;
}
else {
POutput->xx = xx2;
para1 = xx2;
POutput->yy = yy2;
para2 = yy2;
}
POutput->zz = CComplex(1.0,0.0);
}
これで、見た目はうまくいく。px1=(xx1-para1)/(xx2-para1);とおいて、px1.imag()の符号だけ見ればいいかと最初思ったが、px1=0になる可能性があることを忘れていた。一般的には(yy1-para2)/(yy2-para2);の虚部は(xx1-para1)/(xx2-para1);の虚部と同符号だと思うのだが、xx1-xx2がほとんど0の時にはまずいことが起こる。そこで、if(norm(xx1-xx2)>norm(yy1-yy2)){という苦肉の場合分けにした。
(あはら)