diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc index 6bfc9e176644f6543f8d72e59e61e236f5858872..33f8777e4d0a832f2446730c73e928df36bf9406 100644 --- a/tracy/tools/soltracy.cc +++ b/tracy/tools/soltracy.cc @@ -34,8 +34,8 @@ int main(int argc, char *argv[]) { /*set the random value for random generator*/ const long seed = 1121; //the default random seed number - iniranf(seed); //initialize the seed - setrancut(2.0); //default value of the normal cut for the normal distribution + iniranf(seed); //initialize the seed + setrancut(2.0); //default value of the normal cut for the normal distribution // turn on globval.Cavity_on and globval.radiation to get proper synchro radiation damping // IDs accounted too if: wiggler model and symplectic integrator (method = 1) globval.H_exact = false; @@ -67,8 +67,9 @@ int main(int argc, char *argv[]) { if (argc > 1) { - read_script(argv[1], true,CommandNo, UserCommandFlag); - } else { + read_script(argv[1], true, CommandNo, UserCommandFlag); + } + else { fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n"); exit_(1); } @@ -148,11 +149,11 @@ int main(int argc, char *argv[]) { //set RF voltage else if(strcmp(CommandStr,"RFvoltageFlag") == 0) { - printf("\nSetting RF voltage:\n"); - printf(" Old RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav) + fprintf(stdout, "\nSetting RF voltage:\n"); + fprintf(stdout, " Old RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav) / 1e6); set_RFVoltage(globval.cav, UserCommandFlag[i].RFvolt); - printf(" New RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav) + fprintf(stdout, " New RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav) / 1e6); } @@ -293,10 +294,10 @@ int main(int argc, char *argv[]) { qf_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2]; qd_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2]; /* print out the quadrupole strengths before and after the fitting*/ - printf("Before the fitting, the quadrupole field strengths are: \n"); - printf(" %s = %f, %s = %f\n", UserCommandFlag[i].qf, qf_bn, UserCommandFlag[i].qd, qd_bn); - printf("After the fitting, the quadrupole field strengths are: \n"); - printf(" %s = %f, %s = %f\n", UserCommandFlag[i].qf, qf_bn_fit, UserCommandFlag[i].qd, qd_bn_fit); + fprintf(stdout, "Before the fitting, the quadrupole field strengths are: \n"); + fprintf(stdout, " %s = %f, %s = %f\n", UserCommandFlag[i].qf, qf_bn, UserCommandFlag[i].qd, qd_bn); + fprintf(stdout, "After the fitting, the quadrupole field strengths are: \n"); + fprintf(stdout, " %s = %f, %s = %f\n", UserCommandFlag[i].qf, qf_bn_fit, UserCommandFlag[i].qd, qd_bn_fit); /* Compute and get Twiss parameters */ Ring_GetTwiss(chroma = true, 0.0); @@ -333,10 +334,10 @@ int main(int argc, char *argv[]) { qd2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2]; /* print out the quadrupole strengths before and after the fitting*/ printf("Before the fitting, the quadrupole field strengths are: \n"); - printf(" %s = %f, %s = %f, %s = %f, %s = %f\n", UserCommandFlag[i].qf1, qf1_bn, + fprintf(stdout, " %s = %f, %s = %f, %s = %f, %s = %f\n", UserCommandFlag[i].qf1, qf1_bn, UserCommandFlag[i].qf2, qf2_bn, UserCommandFlag[i].qd1, qd1_bn, UserCommandFlag[i].qd2, qd2_bn); printf("After the fitting, the quadrupole field strengths are: \n"); - printf(" %s = %f, %s = %f, %s = %f, %s = %f\n", UserCommandFlag[i].qf1, + fprintf(stdout, " %s = %f, %s = %f, %s = %f, %s = %f\n", UserCommandFlag[i].qf1, qf1_bn_fit, UserCommandFlag[i].qf2, qf2_bn_fit, UserCommandFlag[i].qd1, qd1_bn_fit, UserCommandFlag[i].qd2, qd2_bn_fit); @@ -363,10 +364,10 @@ int main(int argc, char *argv[]) { sxm1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3]; sxm2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3]; /* print out the sextupole strengths before and after the fitting*/ - printf("Before the fitting, the sextupole field strengths are \n"); - printf(" %s = %f, %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn, UserCommandFlag[i].sxm2, sxm2_bn); - printf("After the fitting, the sextupole field strengths are: \n"); - printf(" %s = %f, %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn_fit, UserCommandFlag[i].sxm2, sxm2_bn_fit); + fprintf(stdout, "Before the fitting, the sextupole field strengths are \n"); + fprintf(stdout, " %s = %f, %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn, UserCommandFlag[i].sxm2, sxm2_bn); + fprintf(stdout, "After the fitting, the sextupole field strengths are: \n"); + fprintf(stdout, " %s = %f, %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn_fit, UserCommandFlag[i].sxm2, sxm2_bn_fit); Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */ printglob(); /* print parameter list */ @@ -438,90 +439,119 @@ int main(int argc, char *argv[]) { // Computes FMA else if(strcmp(CommandStr,"FmapFlag") == 0) { - printf("\n begin Fmap calculation for on momentum particles: \n"); + fprintf(stdout, "\n begin Fmap calculation for on momentum particles: \n"); #if MPI_EXEC //initialization for parallel computing - int myid=0, numprocs=0; - int namelen=0; + int myid = 0, numprocs = 0; + int namelen = 0; char processor_name[MPI_MAX_PROCESSOR_NAME]=""; - MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myid); MPI_Get_processor_name(processor_name, &namelen); - printf("\n Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); + fprintf(stdout, "\n Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); + + fprintf(stdout, "\n Begin Fmap calculation for on momentum particles: \n"); - printf("\n begin Fmap calculation for on momentum particles: \n"); - //Each core or process, characterized by myid, will track different region of fmap + //Each core or process, characterized by myid, will track different region of fmap fmap_p(UserCommandFlag[i]._FmapFlag_fmap_file, - UserCommandFlag[i]._FmapFlag_nxpoint, - UserCommandFlag[i]._FmapFlag_nypoint, - UserCommandFlag[i]._FmapFlag_nturn, - UserCommandFlag[i]._FmapFlag_xmax, - UserCommandFlag[i]._FmapFlag_ymax, - UserCommandFlag[i]._FmapFlag_delta, - UserCommandFlag[i]._FmapFlag_diffusion, - UserCommandFlag[i]._FmapFlag_printloss, - numprocs,myid); + UserCommandFlag[i]._FmapFlag_nxpoint, + UserCommandFlag[i]._FmapFlag_nypoint, + UserCommandFlag[i]._FmapFlag_nturn, + UserCommandFlag[i]._FmapFlag_xmax, + UserCommandFlag[i]._FmapFlag_ymax, + UserCommandFlag[i]._FmapFlag_delta, + UserCommandFlag[i]._FmapFlag_diffusion, + UserCommandFlag[i]._FmapFlag_printloss, + numprocs, myid); //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed. MPI_Barrier(MPI_COMM_WORLD); //Collecting data from all files generated by cores, then write to fmap_p.out - if(myid==0) - { - // system("cp 0fmap.out fmap_p.out"); + if(myid == 0){ - FILE *outf; - if ((outf = fopen(UserCommandFlag[i]._FmapFlag_fmap_file, "w")) == NULL) - { - fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file); - exit_(1); - } - - FILE *fp; - char buffer[81]; - char FmapFile[50]; - for(int j=0; j<numprocs; j++) - { - FmapFile[0]='\0'; - sprintf(FmapFile,"%d",j); - strcat(FmapFile,UserCommandFlag[i]._FmapFlag_fmap_file); - - if((fp = fopen(FmapFile, "r")) == NULL) - { - fprintf(stdout, "%s: error while opening file.\n", FmapFile); - exit_(1); - } - - while(fgets(buffer, 80, fp) != NULL) - { - fputs(buffer, outf); - buffer[0]='\0'; - } - fclose(fp); - } - fclose(outf); - } + FILE *outf; + + // open frequency map file for final results + if ((outf = fopen(UserCommandFlag[i]._FmapFlag_fmap_file, "w")) == NULL){ + fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file); + exit_(1); + } + + // open frequency map loss file for final results + if (UserCommandFlag[i]._FmapFlag_printloss){ + FILE *outf_loss; + char FmapLossFile[max_str+5]=" "; + strcpy(FmapLossFile, UserCommandFlag[i]._FmapFlag_fmap_file); + strcat(FmapLossFile, ".loss"); + if ((outf_loss = fopen(FmapLossFile, "w")) == NULL){ + fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file); + exit_(1); + } + } + + FILE *fp, *fp_loss; + char buffer[81]; + char FmapFile[max_str]; + + for(int j = 0; j < numprocs; j++){ + FmapFile[0]='\0'; + sprintf(FmapFile,"%d",j); + strcat(FmapFile, UserCommandFlag[i]._FmapFlag_fmap_file); + + if((fp = fopen(FmapFile, "r")) == NULL) { + fprintf(stdout, "%s: error while opening file.\n", FmapFile); + exit_(1); + while(fgets(buffer, 80, fp) != NULL){ + fputs(buffer, outf); + buffer[0]='\0'; + } + fclose(fp); + } + + if (UserCommandFlag[i]._FmapFlag_printloss){ + strcpy(FmapLossFile,FmapFile); + strcat(FmapLossFile,".loss"); + + if((fp_loss = fopen(FmapFile, "r")) == NULL) { + fprintf(stdout, "%s: error while opening file.\n", FmapFile); + exit_(1); + while(fgets(buffer, 80, fp_loss) != NULL){ + fputs(buffer, outf_loss); + buffer[0]='\0'; + } + } + fclose(fp_loss); + } + } + fclose(outf); + + if (UserCommandFlag[i]._FmapFlag_printloss){ + fclose(outf_loss); + } + } MPI_Barrier(MPI_COMM_WORLD); printf("Process %d finished for on momentum fmap.\n", myid); - #else - fmap(UserCommandFlag[i]._FmapFlag_fmap_file, - UserCommandFlag[i]._FmapFlag_nxpoint, - UserCommandFlag[i]._FmapFlag_nypoint, - UserCommandFlag[i]._FmapFlag_nturn, - UserCommandFlag[i]._FmapFlag_xmax, - UserCommandFlag[i]._FmapFlag_ymax, - UserCommandFlag[i]._FmapFlag_delta, - UserCommandFlag[i]._FmapFlag_diffusion, - UserCommandFlag[i]._FmapFlag_printloss); -#endif + #else // Single processor + fmap(UserCommandFlag[i]._FmapFlag_fmap_file, + UserCommandFlag[i]._FmapFlag_nxpoint, + UserCommandFlag[i]._FmapFlag_nypoint, + UserCommandFlag[i]._FmapFlag_nturn, + UserCommandFlag[i]._FmapFlag_xmax, + UserCommandFlag[i]._FmapFlag_ymax, + UserCommandFlag[i]._FmapFlag_delta, + UserCommandFlag[i]._FmapFlag_diffusion, + UserCommandFlag[i]._FmapFlag_printloss); +#endif // MPI } // Compute FMA dp else if(strcmp(CommandStr,"FmapdpFlag") == 0) { - printf("\n begin Fmap calculation for off momentum particles: \n"); + fprintf(stdout, "\n begin Fmap calculation for off momentum particles: \n"); #if MPI_EXEC @@ -534,70 +564,94 @@ int main(int argc, char *argv[]) { MPI_Get_processor_name(processor_name, &namelen); printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); - printf("\n begin Fmap calculation for off momentum particles: \n"); + fprintf(stdout, "\n begin Fmap calculation for off momentum particles: \n"); fmapdp_p(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, - UserCommandFlag[i]._FmapdpFlag_nxpoint, - UserCommandFlag[i]._FmapdpFlag_nepoint, - UserCommandFlag[i]._FmapdpFlag_nturn, - UserCommandFlag[i]._FmapdpFlag_xmax, - UserCommandFlag[i]._FmapdpFlag_emax, - UserCommandFlag[i]._FmapdpFlag_z, - UserCommandFlag[i]._FmapdpFlag_diffusion, - UserCommandFlag[i]._FmapFlag_printloss, - numprocs,myid); + UserCommandFlag[i]._FmapdpFlag_nxpoint, + UserCommandFlag[i]._FmapdpFlag_nepoint, + UserCommandFlag[i]._FmapdpFlag_nturn, + UserCommandFlag[i]._FmapdpFlag_xmax, + UserCommandFlag[i]._FmapdpFlag_emax, + UserCommandFlag[i]._FmapdpFlag_z, + UserCommandFlag[i]._FmapdpFlag_diffusion, + UserCommandFlag[i]._FmapspFlag_printloss, numprocs, myid); //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed. MPI_Barrier(MPI_COMM_WORLD); //Collecting data from all files generated by cores, then write to fmap_p.out - if(myid==0) - { - // system("cp 0fmap.out fmap_p.out"); - - FILE *outf; - if ((outf = fopen(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, "w")) == NULL) - { - fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmapdp_file); - exit_(1); - } - - FILE *fp; - char buffer[81]; - char FmapFile[50]; - for(int j=0; j<numprocs; j++) - { - FmapFile[0]='\0'; - sprintf(FmapFile,"%d",j); - strcat(FmapFile,UserCommandFlag[i]._FmapdpFlag_fmapdp_file); - - if((fp = fopen(FmapFile, "r")) == NULL) - { - fprintf(stdout, "%s: error while opening file.\n", FmapFile); - exit_(1); - } - - while(fgets(buffer, 80, fp) != NULL) - { - fputs(buffer, outf); - buffer[0]='\0'; - } - fclose(fp); - } - fclose(outf); - } - MPI_Barrier(MPI_COMM_WORLD); - printf("Process %d finished off momentum fmap.\n", myid); - - #else - fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, - UserCommandFlag[i]._FmapdpFlag_nxpoint, - UserCommandFlag[i]._FmapdpFlag_nepoint, - UserCommandFlag[i]._FmapdpFlag_nturn, - UserCommandFlag[i]._FmapdpFlag_xmax, - UserCommandFlag[i]._FmapdpFlag_emax, - UserCommandFlag[i]._FmapdpFlag_z, - UserCommandFlag[i]._FmapdpFlag_diffusion, - UserCommandFlag[i]._FmapdpFlag_printloss); + if(myid == 0){ + FILE *outf; + if ((outf = fopen(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, "w")) == NULL){ + fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmapdp_file); + exit_(1); + } + + // open frequency map loss file for final results + if (UserCommandFlag[i]._FmapdpFlag_printloss){ + FILE *outf_loss; + char FmapdpLossFile[max_str+5]=" "; + strcpy(FmapdpLossFile, UserCommandFlag[i]._FmapdpFlag_fmapdp_file); + strcat(FmapdpLossFile, ".loss"); + if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL){ + fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmap_file); + exit_(1); + } + } + + FILE *fp, *fp_loss; + char buffer[81]; + char FmapFile[50]; + for(int j = 0; j < numprocs; j++){ + FmapFile[0]='\0'; + sprintf(FmapFile,"%d",j); + strcat(FmapFile,UserCommandFlag[i]._FmapdpFlag_fmapdp_file); + + if((fp = fopen(FmapFile, "r")) == NULL){ + fprintf(stdout, "%s: error while opening file.\n", FmapFile); + exit_(1); + } + + while(fgets(buffer, 80, fp) != NULL){ + fputs(buffer, outf); + buffer[0]='\0'; + } + fclose(fp); + + if (UserCommandFlag[i]._FmapdpFlag_printloss){ + strcpy(FmapdpLossFile, FmapdpFile); + strcat(FmapdpLossFile,".loss"); + + if((fp_loss = fopen(FmapdpFile, "r")) == NULL) { + fprintf(stdout, "%s: error while opening file.\n", FmapdpFile); + exit_(1); + while(fgets(buffer, 80, fp_loss) != NULL){ + fputs(buffer, outf_loss); + buffer[0]='\0'; + } + } + fclose(fp_loss); + } // if printloss + } // for loop + fclose(outf); + + if (UserCommandFlag[i]._FmapFlag_printloss){ + fclose(outf_loss); + } + } // if myid + + MPI_Barrier(MPI_COMM_WORLD); + fprintf(stdout, "Process %d finished off momentum fmap.\n", myid); + + #else //single processor + fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, + UserCommandFlag[i]._FmapdpFlag_nxpoint, + UserCommandFlag[i]._FmapdpFlag_nepoint, + UserCommandFlag[i]._FmapdpFlag_nturn, + UserCommandFlag[i]._FmapdpFlag_xmax, + UserCommandFlag[i]._FmapdpFlag_emax, + UserCommandFlag[i]._FmapdpFlag_z, + UserCommandFlag[i]._FmapdpFlag_diffusion, + UserCommandFlag[i]._FmapdpFlag_printloss); #endif } @@ -626,7 +680,7 @@ int main(int argc, char *argv[]) { }; /* calculate momentum acceptance*/ - printf("\n Calculate momentum acceptance: \n"); + fprintf(stdout, "\n Calculate momentum acceptance: \n"); #if MPI_EXEC @@ -640,126 +694,112 @@ int main(int argc, char *argv[]) { MPI_Get_processor_name(processor_name, &namelen); printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); - printf("\n Calculate momentum acceptance: \n"); + fprintf(stdout, "\n Calculate momentum acceptance: \n"); MomentumAcceptance_p(UserCommandFlag[i]._MomentumAccFlag_momacc_file, - UserCommandFlag[i]._MomentumAccFlag_istart, - UserCommandFlag[i]._MomentumAccFlag_istop, - UserCommandFlag[i]._MomentumAccFlag_deltaminp, - UserCommandFlag[i]._MomentumAccFlag_deltamaxp, - UserCommandFlag[i]._MomentumAccFlag_nstepp, - UserCommandFlag[i]._MomentumAccFlag_deltaminn, - UserCommandFlag[i]._MomentumAccFlag_deltamaxn, - UserCommandFlag[i]._MomentumAccFlag_nstepn, - UserCommandFlag[i]._MomentumAccFlag_nturn, - UserCommandFlag[i]._MomentumAccFlag_zmax, - numprocs,myid); + UserCommandFlag[i]._MomentumAccFlag_istart, + UserCommandFlag[i]._MomentumAccFlag_istop, + UserCommandFlag[i]._MomentumAccFlag_deltaminp, + UserCommandFlag[i]._MomentumAccFlag_deltamaxp, + UserCommandFlag[i]._MomentumAccFlag_nstepp, + UserCommandFlag[i]._MomentumAccFlag_deltaminn, + UserCommandFlag[i]._MomentumAccFlag_deltamaxn, + UserCommandFlag[i]._MomentumAccFlag_nstepn, + UserCommandFlag[i]._MomentumAccFlag_nturn, + UserCommandFlag[i]._MomentumAccFlag_zmax, + numprocs,myid); //merge the files //Synchronize all cores, till finish cal of different region, //then files from the cores will be processed. MPI_Barrier(MPI_COMM_WORLD); //Collect data from all files generated by different cores, then write to user defined file - if(myid==0) - { - //merge the phase.out from different cpus - FILE *outf1; //phase.out, for debugging - FILE *fp1; - char buffer1[81]; - char PhaseFile[50]; - if ((outf1 = fopen("phase_p.out", "w")) == NULL) - { - fprintf(stdout, "psoltracy: error while opening file %s\n", "phase_p.out"); - exit_(1); - } - - for(int j=0; j<numprocs; j++) - { - PhaseFile[0]='\0'; - sprintf(PhaseFile,"%d",j); - strcat(PhaseFile,"phase.out"); - - if((fp1 = fopen(PhaseFile, "r")) == NULL){ - fprintf(stdout, "%s: error while opening file.\n", PhaseFile); - exit_(1); - } - - while(fgets(buffer1, 80, fp1) != NULL) - { - fputs(buffer1, outf1); - buffer1[0]='\0'; - } - fclose(fp1); - } - - //collect data for the momentum acceptance - FILE *outf2; //momentum acceptance - FILE *fp2; - char buffer2[81]; - char MAFile[50]; - if ((outf2 = fopen(UserCommandFlag[i]._MomentumAccFlag_momacc_file, "w")) == NULL) - { - fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._MomentumAccFlag_momacc_file); - exit_(1); - } - - //Collect data in Positive region - for(int j=0; j<numprocs; j++) - { - MAFile[0]='\0'; - sprintf(MAFile,"%d",j); - strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file); - - if((fp2 = fopen(MAFile, "r")) == NULL) - { - fprintf(stdout, "%s: error while opening file.\n", MAFile); - exit_(1); - } + if(myid==0){ + //merge the phase.out from different cpus + FILE *outf1; //phase.out, for debugging + FILE *fp1; + char buffer1[81]; + char PhaseFile[50]; + if ((outf1 = fopen("phase_p.out", "w")) == NULL){ + fprintf(stdout, "psoltracy: error while opening file %s\n", "phase_p.out"); + exit_(1); + } + + for(int j=0; j<numprocs; j++){ + PhaseFile[0]='\0'; + sprintf(PhaseFile,"%d",j); + strcat(PhaseFile,"phase.out"); + if((fp1 = fopen(PhaseFile, "r")) == NULL){ + fprintf(stdout, "%s: error while opening file.\n", PhaseFile); + exit_(1); + } + + while(fgets(buffer1, 80, fp1) != NULL){ + fputs(buffer1, outf1); + buffer1[0]='\0'; + } + fclose(fp1); + } // for loop + + //collect data for the momentum acceptance + FILE *outf2; //momentum acceptance + FILE *fp2; + char buffer2[81]; + char MAFile[50]; + if ((outf2 = fopen(UserCommandFlag[i]._MomentumAccFlag_momacc_file, "w")) == NULL){ + fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._MomentumAccFlag_momacc_file); + exit_(1); + } + + //Collect data in Positive region + for(int j=0; j<numprocs; j++){ + MAFile[0]='\0'; + sprintf(MAFile,"%d",j); + strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file); + + if((fp2 = fopen(MAFile, "r")) == NULL){ + fprintf(stdout, "%s: error while opening file.\n", MAFile); + exit_(1); + } - while(fgets(buffer2, 80, fp2) != NULL) - { - if(strstr(buffer2,"Negative")!=0) break; - - fputs(buffer2, outf2); - buffer2[0]='\0'; //reset buffer to NULL - } - buffer2[0]='\0'; - fclose(fp2); - } - - fprintf(outf2,"\n"); // A void line - - //Collect data in Nagative region - for(int j=0; j<numprocs; j++) - { - MAFile[0]='\0'; - sprintf(MAFile,"%d",j); - strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file); - - if((fp2 = fopen(MAFile, "r")) == NULL) - { - fprintf(stdout, "%s: error while opening file.\n", MAFile); - exit_(1); - } - - bool found=false; - while(fgets(buffer2, 80, fp2) != NULL) - { - if( (strstr(buffer2,"Negative")==0) && (!found)) { buffer2[0]='\0'; continue;} - if( (strstr(buffer2,"Negative")!=0) ) { found=true; buffer2[0]='\0'; continue;} - - fputs(buffer2, outf2); - buffer2[0]='\0'; - } - buffer2[0]='\0'; - fclose(fp2); - } - fclose(outf2); - - } + while(fgets(buffer2, 80, fp2) != NULL){ + if(strstr(buffer2,"Negative")!=0) break; + fputs(buffer2, outf2); + buffer2[0]='\0'; //reset buffer to NULL + } + buffer2[0]='\0'; + fclose(fp2); + } + + fprintf(outf2,"\n"); // A void line + + //Collect data in Nagative region + for(int j=0; j<numprocs; j++){ + MAFile[0]='\0'; + sprintf(MAFile,"%d",j); + strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file); + + if((fp2 = fopen(MAFile, "r")) == NULL){ + fprintf(stdout, "%s: error while opening file.\n", MAFile); + exit_(1); + } + + bool found=false; + while(fgets(buffer2, 80, fp2) != NULL){ + if((strstr(buffer2,"Negative")==0) && (!found)) { buffer2[0]='\0'; continue;} + if((strstr(buffer2,"Negative")!=0) ) { found=true; buffer2[0]='\0'; continue;} + + fputs(buffer2, outf2); + buffer2[0]='\0'; + } // while + buffer2[0]='\0'; + fclose(fp2); + } // for + fclose(outf2); + } MPI_Barrier(MPI_COMM_WORLD); printf("process %d finished momentum acceptance.\n", myid); -#else +#else // single processor coputation MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file, UserCommandFlag[i]._MomentumAccFlag_istart, UserCommandFlag[i]._MomentumAccFlag_istop, @@ -780,11 +820,10 @@ int main(int argc, char *argv[]) { // induced amplitude else if(strcmp(CommandStr,"InducedAmplitudeFlag") == 0) { - printf("\n Calculate induced amplitude: \n"); - InducedAmplitude(193L); + fprintf(stdout, "\n Calculate induced amplitude: \n"); + InducedAmplitude(193L); } - else if(strcmp(CommandStr,"EtaFlag") == 0) { // compute cod and Twiss parameters for different energy offsets for (int ii = 0; ii <= 40; ii++) { @@ -829,14 +868,14 @@ int main(int argc, char *argv[]) { Phase(UserCommandFlag[i]._Phase_phase_file, UserCommandFlag[i]._Phase_X, UserCommandFlag[i]._Phase_Px, - UserCommandFlag[i]._Phase_Y, - UserCommandFlag[i]._Phase_Py, - UserCommandFlag[i]._Phase_delta, - UserCommandFlag[i]._Phase_ctau, - UserCommandFlag[i]._Phase_nturn); + UserCommandFlag[i]._Phase_Y, + UserCommandFlag[i]._Phase_Py, + UserCommandFlag[i]._Phase_delta, + UserCommandFlag[i]._Phase_ctau, + UserCommandFlag[i]._Phase_nturn); printf("6D phase space tracking: \n the simulation time for phase space in tracy 3 is \n"); stop = stampstop(start); - + /* restore the initial values*/ globval.Cavity_on = cavityflag; globval.radiation = radiationflag; @@ -966,7 +1005,7 @@ Based on parts of functions get_param( ) & ID_corr(), etc in nsls-ii_lib.cc cout << "Number of quadrupole families used for ID correction: " << N_Fam << endl; cout << "Quadrupoles used for ID correction: " << endl; for (int k = 0; k < N_Fam; k++) - printf("%d\n",Q_Fam[k]); + fprintf(stdout, "%d\n",Q_Fam[k]); } //initialization diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc index ef7d4619770ebebc50c15bb01b5f3b3caf7fc5f4..dc50e7ba9c81599a4b2cf26d83d3c24776c318d1 100644 --- a/tracy/tracy/src/read_script.cc +++ b/tracy/tracy/src/read_script.cc @@ -35,30 +35,30 @@ char girder_file[max_str]; - // multipole files; for soleil lattice +// multipole files; for soleil lattice char multipole_file[max_str]; - // multipole file for soleil lattice +// multipole file for soleil lattice char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str]; - // files to set sources of coupling; for SOLEIL lattice +// files to set sources of coupling; for SOLEIL lattice char virtualskewquad_file[max_str]; //errors //char fe_file[max_str]; //the same as multipole_file[max_str]???? /* COD correction */ -int nwh = 56, nwv = 56; //number of singular values for SVD correction +int nwh = 57, nwv = 57; //number of singular values for SVD correction /* ID compensation */ -char IDCq_name[max_str][11]; +char IDCq_name[max_str][11]; //multipole field error correction (SOLEIL based) char hcorr_file[max_str], vcorr_file[max_str]; //files with the status of hcorr/vcorr status, //to choose which correctors are used for orbit correction //default name of special elements in lattice -char hcorr_name[max_str] = "", vcorr_name[max_str] = ""; -char skew_quad_name[max_str] = "", bpm_name[max_str] = ""; -char gs_name[max_str] = "", ge_name[max_str] = ""; +char hcorr_name[max_str] = "", vcorr_name[max_str] = ""; +char skew_quad_name[max_str] = "", bpm_name[max_str] = ""; +char gs_name[max_str] = "", ge_name[max_str] = ""; #define OLD_LATTICE @@ -226,78 +226,76 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom //print close orbit(COD) flag else if (strcmp("PrintCODFlag", name) == 0){ sscanf(line, "%*s %s",nextpara); - - if(strcmp(nextpara,"voidpara")!=0) - sscanf(line, "%*s %s", UserCommandFlag[CommNo]._PrintCOD_cod_file); - - strcpy(UserCommandFlag[CommNo].CommandStr,name); + if(strcmp(nextpara,"voidpara")!=0) + sscanf(line, "%*s %s", UserCommandFlag[CommNo]._PrintCOD_cod_file); + strcpy(UserCommandFlag[CommNo].CommandStr,name); } //print the cooridinates using tracking at each element else if (strcmp("PrintTrackFlag", name) == 0){ - sscanf(line, "%*s %s",nextpara); + sscanf(line, "%*s %s",nextpara); - if(strcmp(nextpara,"voidpara")!=0) - sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld", - UserCommandFlag[CommNo]._PrintTrack_track_file, - &(UserCommandFlag[CommNo]._PrintTrack_x), - &(UserCommandFlag[CommNo]._PrintTrack_px), - &(UserCommandFlag[CommNo]._PrintTrack_y), - &(UserCommandFlag[CommNo]._PrintTrack_py), - &(UserCommandFlag[CommNo]._PrintTrack_delta), - &(UserCommandFlag[CommNo]._PrintTrack_ctau), - &(UserCommandFlag[CommNo]._PrintTrack_nmax)); - strcpy(UserCommandFlag[CommNo].CommandStr,name); - } - //print close orbit(COD) flag + if(strcmp(nextpara,"voidpara")!=0) + sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld", + UserCommandFlag[CommNo]._PrintTrack_track_file, + &(UserCommandFlag[CommNo]._PrintTrack_x), + &(UserCommandFlag[CommNo]._PrintTrack_px), + &(UserCommandFlag[CommNo]._PrintTrack_y), + &(UserCommandFlag[CommNo]._PrintTrack_py), + &(UserCommandFlag[CommNo]._PrintTrack_delta), + &(UserCommandFlag[CommNo]._PrintTrack_ctau), + &(UserCommandFlag[CommNo]._PrintTrack_nmax)); + strcpy(UserCommandFlag[CommNo].CommandStr,name); + } + + //print close orbit(COD) flag else if (strcmp("PrintGirderFlag", name) == 0){ - sscanf(line, "%*s %s", girder_file); - strcpy(UserCommandFlag[CommNo].CommandStr,name); + sscanf(line, "%*s %s", girder_file); + strcpy(UserCommandFlag[CommNo].CommandStr,name); } else if (strcmp("TuneTracFlag", name) == 0){ - strcpy(UserCommandFlag[CommNo].CommandStr,name); + strcpy(UserCommandFlag[CommNo].CommandStr,name); } else if (strcmp("ChromTracFlag", name) == 0){ - strcpy(UserCommandFlag[CommNo].CommandStr,name); + strcpy(UserCommandFlag[CommNo].CommandStr,name); } - // FMA - else if (strcmp("FmapFlag", name) == 0){ - strcpy(dummy, ""); - strcpy(dummy2,""); - sscanf(line, "%*s %s",nextpara); + //FMA + else if (strcmp("FmapFlag", name) == 0){ + strcpy(dummy , ""); + strcpy(dummy2, ""); + sscanf(line, "%*s %s",nextpara); - //if no definition of loss flag in the lattice file, then "dummy2" is empty string. - if(strcmp(nextpara,"voidpara")!=0) - sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %s", - UserCommandFlag[CommNo]._FmapFlag_fmap_file, - &(UserCommandFlag[CommNo]._FmapFlag_nxpoint), - &(UserCommandFlag[CommNo]._FmapFlag_nypoint), - &(UserCommandFlag[CommNo]._FmapFlag_nturn), - &(UserCommandFlag[CommNo]._FmapFlag_xmax), - &(UserCommandFlag[CommNo]._FmapFlag_ymax), - &(UserCommandFlag[CommNo]._FmapFlag_delta), - dummy,dummy2); - - if(strcmp(dummy, "true") == 0) - UserCommandFlag[CommNo]._FmapFlag_diffusion = true; - else if(strcmp(dummy, "false") == 0) - UserCommandFlag[CommNo]._FmapFlag_diffusion = false; - - if(strcmp(dummy2, "true") == 0) - UserCommandFlag[CommNo]._FmapFlag_printloss = true; - else if(strcmp(dummy2, "false") == 0) - UserCommandFlag[CommNo]._FmapFlag_printloss = false; - else - UserCommandFlag[CommNo]._FmapFlag_printloss = false; - - - //cout << "debug: " << line << endl; - //cout << "dummy2 = " << dummy2 << " lossflag = " << UserCommandFlag[CommNo]._FmapFlag_printloss << endl; - - // FmapFlag = true; - strcpy(UserCommandFlag[CommNo].CommandStr,name); - } + //if no definition of loss flag in the lattice file, then "dummy2" is empty string. + if(strcmp(nextpara,"voidpara")!= 0) + sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %s", + UserCommandFlag[CommNo]._FmapFlag_fmap_file, + &(UserCommandFlag[CommNo]._FmapFlag_nxpoint), + &(UserCommandFlag[CommNo]._FmapFlag_nypoint), + &(UserCommandFlag[CommNo]._FmapFlag_nturn), + &(UserCommandFlag[CommNo]._FmapFlag_xmax), + &(UserCommandFlag[CommNo]._FmapFlag_ymax), + &(UserCommandFlag[CommNo]._FmapFlag_delta), + dummy,dummy2); + + if(strcmp(dummy, "true") == 0) + UserCommandFlag[CommNo]._FmapFlag_diffusion = true; + else if(strcmp(dummy, "false") == 0) + UserCommandFlag[CommNo]._FmapFlag_diffusion = false; + + if(strcmp(dummy2, "true") == 0) + UserCommandFlag[CommNo]._FmapFlag_printloss = true; + else if(strcmp(dummy2, "false") == 0) + UserCommandFlag[CommNo]._FmapFlag_printloss = false; + else + UserCommandFlag[CommNo]._FmapFlag_printloss = false; + + //cout << "debug: " << line << endl; + //cout << "dummy2 = " << dummy2 << " lossflag = " << UserCommandFlag[CommNo]._FmapFlag_printloss << endl; + + // FmapFlag = true; + strcpy(UserCommandFlag[CommNo].CommandStr,name); + } // FMA dp else if (strcmp("FmapdpFlag", name) == 0){ strcpy(dummy, ""); @@ -357,7 +355,7 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom sscanf(line, "%*s %s %ld %ld %lf", UserCommandFlag[CommNo]._EnergyTuneShift_nudp_file, &(UserCommandFlag[CommNo]._EnergyTuneShift_npoint), - &(UserCommandFlag[CommNo]._EnergyTuneShift_nturn), + &(UserCommandFlag[CommNo]._EnergyTuneShift_nturn), &(UserCommandFlag[CommNo]._EnergyTuneShift_deltamax)); strcpy(UserCommandFlag[CommNo].CommandStr,name); @@ -386,13 +384,13 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom UserCommandFlag[CommNo]._MomentumAccFlag_TrackDim, &(UserCommandFlag[CommNo]._MomentumAccFlag_istart), &(UserCommandFlag[CommNo]._MomentumAccFlag_istop), - &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminp), + &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminp), &(UserCommandFlag[CommNo]._MomentumAccFlag_deltamaxp), &(UserCommandFlag[CommNo]._MomentumAccFlag_nstepp), - &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminn), + &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminn), &(UserCommandFlag[CommNo]._MomentumAccFlag_deltamaxn), &(UserCommandFlag[CommNo]._MomentumAccFlag_nstepn), - &(UserCommandFlag[CommNo]._MomentumAccFlag_nturn), + &(UserCommandFlag[CommNo]._MomentumAccFlag_nturn), &(UserCommandFlag[CommNo]._MomentumAccFlag_zmax)); } strcpy(UserCommandFlag[CommNo].CommandStr,name); @@ -406,15 +404,15 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom } // fit for 2 families,specific for the soleil lattice in which the quadrupole is cut into 2 parts else if (strcmp("FitTune4Flag", name) == 0){ - sscanf(line, "%*s %s %s %s %s %lf %lf",UserCommandFlag[CommNo].qf1,UserCommandFlag[CommNo].qf2, - UserCommandFlag[CommNo].qd1,UserCommandFlag[CommNo].qd2, - &(UserCommandFlag[CommNo].targetnux),&(UserCommandFlag[CommNo].targetnuz)); - strcpy(UserCommandFlag[CommNo].CommandStr,name); + sscanf(line, "%*s %s %s %s %s %lf %lf",UserCommandFlag[CommNo].qf1,UserCommandFlag[CommNo].qf2, + UserCommandFlag[CommNo].qd1,UserCommandFlag[CommNo].qd2, + &(UserCommandFlag[CommNo].targetnux),&(UserCommandFlag[CommNo].targetnuz)); + strcpy(UserCommandFlag[CommNo].CommandStr,name); } else if (strcmp("FitChromFlag", name) == 0){ - sscanf(line, "%*s %s %s %lf %lf",UserCommandFlag[CommNo].sxm1,UserCommandFlag[CommNo].sxm2, - &(UserCommandFlag[CommNo].targetksix),&(UserCommandFlag[CommNo].targetksiz)); - strcpy(UserCommandFlag[CommNo].CommandStr,name); + sscanf(line, "%*s %s %s %lf %lf",UserCommandFlag[CommNo].sxm1,UserCommandFlag[CommNo].sxm2, + &(UserCommandFlag[CommNo].targetksix),&(UserCommandFlag[CommNo].targetksiz)); + strcpy(UserCommandFlag[CommNo].CommandStr,name); } // else if (strcmp("GirderErrorFlag", name) == 0){ // strcpy(UserCommand[CommandNo-1],name); @@ -432,7 +430,7 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom // // } else if (strcmp("InducedAmplitudeFlag", name) == 0){ - strcpy(UserCommandFlag[CommNo].CommandStr,name); + strcpy(UserCommandFlag[CommNo].CommandStr,name); } // else if (strcmp("CodeComparaisonFlag", name) == 0){ // strcpy(UserCommand[CommandNo-1],name); @@ -527,27 +525,27 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom else if (strcmp("IDCorrFlag", name) == 0) strcpy(UserCommandFlag[CommNo].CommandStr,name); else if (strcmp("scl_dbetax", name) == 0)//scaling of dbetax, for ID correction - sscanf(line, "%*s %lf", &scl_dbetax); + sscanf(line, "%*s %lf", &scl_dbetax); else if (strcmp("scl_dbetay", name) == 0)//scaling of dbetay, for ID correction - sscanf(line, "%*s %lf", &scl_dbetay); + sscanf(line, "%*s %lf", &scl_dbetay); else if (strcmp("scl_dnux", name) == 0)//scaling of dnux, for ID correction - sscanf(line, "%*s %lf", &scl_dnux); + sscanf(line, "%*s %lf", &scl_dnux); else if (strcmp("scl_dnuy", name) == 0)//scaling of dnuy, for ID correction - sscanf(line, "%*s %lf", &scl_dnuy); + sscanf(line, "%*s %lf", &scl_dnuy); else if (strcmp("scl_nux", name) == 0) //scaling of nux, for ID correction - sscanf(line, "%*s %lf", &scl_nux); + sscanf(line, "%*s %lf", &scl_nux); else if (strcmp("scl_nuy", name) == 0) //scaling of nuy, for ID correction - sscanf(line, "%*s %lf", &scl_nuy); + sscanf(line, "%*s %lf", &scl_nuy); else if (strcmp("ID_step", name) == 0)//ID steps, for ID correction - sscanf(line, "%*s %lf", &ID_step); + sscanf(line, "%*s %lf", &ID_step); else if (strcmp("N_calls", name) == 0) // ID correction parameters - sscanf(line, "%*s %d", &N_calls); + sscanf(line, "%*s %d", &N_calls); else if (strcmp("N_steps", name) == 0) - sscanf(line, "%*s %d", &N_steps); + sscanf(line, "%*s %d", &N_steps); else if (strcmp("N_Fam", name) == 0) - sscanf(line, "%*s %d", &N_Fam); + sscanf(line, "%*s %d", &N_Fam); else if (strcmp("IDCquads", name) == 0) { - sscanf(line, "%*s %s %s %s %s %s %s %s %s %s %s %s", + sscanf(line, "%*s %s %s %s %s %s %s %s %s %s %s %s", IDCq_name[0], IDCq_name[1], IDCq_name[2], IDCq_name[3], IDCq_name[4], IDCq_name[5], IDCq_name[6], IDCq_name[7], IDCq_name[8], IDCq_name[9], IDCq_name[10]); diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc index 1f6f106836dc383bb6dc2c842881a67523de4a62..705fd2fed0fab4ddd284e6ab77310759dcfbb2ba 100644 --- a/tracy/tracy/src/soleillib.cc +++ b/tracy/tracy/src/soleillib.cc @@ -223,7 +223,7 @@ Vector2 H; at position pos and for a energy offset dP For a givien delta - H = gamma xcod² + 2*alpha*xcod*xcod' + beta*xcod'*xcod' + H = gamma xcod² + 2*alpha*xcod*xcod' + beta*xcod'*xcod' Input: none @@ -875,7 +875,7 @@ void fmap(const char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, do FILE *outf_loss; //file with the information of the lost particle long i = 0L, j = 0L; double Tab[DIM][NTURN], Tab0[DIM][NTURN]; - double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx, dfy; + double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx = 0.0, dfy = 0.0; double x = 0.0, xp = 0.0, y = 0.0, zp = 0.0; double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0; const double ctau = 0.0; @@ -892,12 +892,11 @@ void fmap(const char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, do long lastpos = 1; ss_vect<double> x1; - - //if(loss){ - char FmapLossFile[S_SIZE+5]=" "; - strcpy(FmapLossFile,FmapFile); - strcat(FmapLossFile,".loss"); - //} + if(printloss){ + char FmapLossFile[S_SIZE+5]=" "; + strcpy(FmapLossFile,FmapFile); + strcat(FmapLossFile,".loss"); + } /* Get time and date */ time_t aclock; @@ -1081,15 +1080,15 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax char FmapFile[max_str]; sprintf(FmapFile,"%d",myid); strcat(FmapFile,FmapFile_p); - fprintf(stdout, "%s\n",FmapFile); + if (trace) fprintf(stdout, "%s\n",FmapFile); //variables of the returned the tracked particle - long lastn = 0; - long lastpos = 1; + long lastn = 0L; + long lastpos = 1L; ss_vect<double> x1; char FmapLossFile[S_SIZE+5]=" "; - strcpy(FmapLossFile,FmapFile); + strcpy(FmapLossFile, FmapFile); strcat(FmapLossFile,".loss"); @@ -1484,7 +1483,7 @@ void fmapdp_p(const char *FmapdpFile_p, long Nbx, long Nbe, long Nbtour, double long lastpos = 1; ss_vect<double> x1; - char FmapdpLossFile[S_SIZE+5]=" "; + char FmapdpLossFile[max_str+5]=" "; strcpy(FmapdpLossFile,FmapdpFile); strcat(FmapdpLossFile,".loss");