AUTO から,shooting point を得る.
wrtbv8(iap_type *iap, rap_type *rap, doublereal *par, integer *icp, doublereal *rldot, integer *ndxloc, doublereal **ups, doublereal **udotps, doublereal *tm, doublereal *dtm) { ... mtot = ntot % 10000; fprintf(fp8,"%5ld",ibr); fprintf(fp8,"%5ld",mtot); fprintf(fp8,"%5ld",itp); fprintf(fp8,"%5ld",lab); fprintf(fp8,"%5ld",nfpr); fprintf(fp8,"%5ld",isw); fprintf(fp8,"%5ld",ntpl); fprintf(fp8,"%5ld",nar); fprintf(fp8,"%7ld",nrowpr); fprintf(fp8,"%5ld",ntst); fprintf(fp8,"%5ld",ncol); fprintf(fp8,"%5d\n",NPARX); // この下を追記. // LABEL XXX の時のメッシュ点を splXXX.dat を出力する. //======================================================== // uniker : output the data of the orbit // ntst hooting points are got from AUTO // v[0],..., v[ntst-1] // v[ntst] sim v[0] char char_fname_sp[20]; sprintf(char_fname_sp, "spl%03d.dat", mtot); FILE *fsp=fopen(char_fname_sp,"w"); if(fsp == NULL) { printf("# ERROR: Could not open spXXX.dat"); abort(); } double period = par[10]; //tm[j] in [0.0, 1.0] double s_j; // in [0.0, T] double h_j; // = s[j+1] - s[j] for (j = 0; j < ntst; ++j) { s_j = period*tm[j]; h_j = period*(tm[j+1] - tm[j]); fprintf(fsp," %ld", j); fprintf(fsp,"%19.10E", tm[j]); fprintf(fsp,"%19.10E", s_j); fprintf(fsp,"%19.10E", h_j); for (i = 0; i < ndim; ++i) { fprintf(fsp,"%19.10E", ups[j][i]); } fprintf(fsp,"\n"); } s_j = period*tm[ntst]; // period*1.0 h_j = period - s_j; // 0.0 fprintf(fsp," %ld", ntst); fprintf(fsp,"%19.10E", tm[ntst]); fprintf(fsp,"%19.10E", s_j); fprintf(fsp,"%19.10E", h_j); for (i = 0; i < ndim; ++i) { fprintf(fsp,"%19.10E", ups[ntst][i]); } fprintf(fsp,"\n"); //======================================================== // ここまで.この後は,元々の処理が続く /* Write the entire solution on unit 8 : */ for (j = 0; j < ntst; ++j) { rn = 1. / ncol; for (i = 0; i < ncol; ++i) { k1 = i * ndim; k2 = (i + 1) * ndim; t = tm[j] + i * rn * dtm[j]; fprintf(fp8," %19.10E",t); ...
なお,ups の中身は,以下のようになっている.
// example // dimension = 2 // 0<= j < ntst // (ncol-1) interpolated points //ups[j][0] // j_0 //ups[j][1] //ups[j][2] // j_1 //ups[j][3] //... //ups[j][(ncol-1)*dim] // j_{ncol-1} //ups[j][(ncol-1)*dim + 1] //ups[j+1][0] //ups[j+1][1] //... //ups[ntst][0] //ups[ntst][1]