とりあえず、後半部分をこうしてみた。
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)){という苦肉の場合分けにした。
(あはら)