From 08783d4bb5157f705dcfc16608c4ad200312c0b1 Mon Sep 17 00:00:00 2001 From: Tobias Date: Thu, 19 Mar 2020 18:08:35 +0100 Subject: [PATCH 01/20] This seems to give a working distlib.. --- CMakeLists.txt | 1 + src/CMakeLists.txt | 1 + src/mad_cmd.c | 2 ++ src/mad_dict.c | 10 ++++++++++ src/mad_extrn_f.h | 1 + src/mad_gcst.c | 16 +++++++++++++++- src/mad_gcst.h | 3 +++ src/madx.h | 3 ++- 8 files changed, 35 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f6922ec6b..1ba4cd1f9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,6 +30,7 @@ include(setupCompilerSpecifics) include(setupArchSpecifics) add_subdirectory(libs/ptc) +add_subdirectory(libs/DISTlib) add_subdirectory(src) add_subdirectory(tools) add_subdirectory(syntax) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3ce37764d..4c447eb34 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -33,6 +33,7 @@ file(GLOB headerfiles *.h) add_library(madx ${src_files}) target_link_libraries(madx ptc) +target_link_libraries(madx DISTlib) if(USE_GC) target_link_libraries(madx gc-lib) diff --git a/src/mad_cmd.c b/src/mad_cmd.c index c96250801..6623a7c91 100644 --- a/src/mad_cmd.c +++ b/src/mad_cmd.c @@ -181,6 +181,7 @@ exec_command(void) current_command = p->clone; // main commands + printf("outtt %s \n", p->cmd_def->module ); if (strcmp(p->cmd_def->module, "control") == 0) control(p); else if (strcmp(p->cmd_def->module, "c6t") == 0) conv_sixtrack(p); else if (strcmp(p->cmd_def->module, "edit") == 0) seq_edit_main(p); @@ -196,6 +197,7 @@ exec_command(void) else if (strcmp(p->cmd_def->module, "survey") == 0) { current_survey = p->clone; pro_survey(p); } else if (strcmp(p->cmd_def->module, "track") == 0) pro_track(p); else if (strcmp(p->cmd_def->module, "twiss") == 0) { current_twiss = p->clone; pro_twiss(); } + else if (strcmp(p->cmd_def->module, "distribution") == 0) pro_distribution(p); else if (strcmp(p->cmd_def->module, "sdds") == 0) { #ifdef _ONLINE diff --git a/src/mad_dict.c b/src/mad_dict.c index ba013f1b5..de9963728 100644 --- a/src/mad_dict.c +++ b/src/mad_dict.c @@ -230,6 +230,16 @@ const char *const_command_def = "sequence = [s, none], " "table = [s, none]; " " " +"distribution: distribution none 0 0 " +"npart = [i, -1], " +"horizontal = [s, zero, zero], " +"vertical = [s, zero, zero], " +"longitudinal = [s, zero, zero], " +"table = [s, disttable, disttable], " +"cutsig_h = [r, {-1}], " +"cutsig_v = [r, {-1}], " +"cutsig_l = [r, {-1}]; " +" " "deselect: control none 0 0 " /* GR: stub exists but not properly implemented */ "flag = [s, none, none], " "range = [s, #s/#e, none], " diff --git a/src/mad_extrn_f.h b/src/mad_extrn_f.h index 9edcf6adc..6c5f175f1 100644 --- a/src/mad_extrn_f.h +++ b/src/mad_extrn_f.h @@ -283,6 +283,7 @@ void trrun_(F_INTEGER switch_, F_INTEGER turns, F_DOUBLE orbit0, F_DOUBLE rt, F_ F_INTEGER last_turn, F_DOUBLE last_pos, F_DOUBLE z, F_DOUBLE dxt, F_DOUBLE dyt, F_DOUBLE last_orbit, F_DOUBLE eigen, F_DOUBLE coords, F_INTEGER e_flag, F_INTEGER code_buf, F_DOUBLE l_buf); +void trbegn_(F_DOUBLE rt, F_DOUBLE eigen); // from twiss.f90 void tmrefe_(F_DOUBLE rf); diff --git a/src/mad_gcst.c b/src/mad_gcst.c index 29f57195d..153512622 100644 --- a/src/mad_gcst.c +++ b/src/mad_gcst.c @@ -113,12 +113,14 @@ const char* const survey_table_cols[] = "name", "keyword", "s", "l", "angle", "x", "y", "z", "theta", "phi", "psi", "globaltilt", "tilt", "slot_id", "assembly_id", "mech_sep", - /*== dealt with the new property v_pos as for mech_sep */ + /*== dealt with the new property v_pos asc for mech_sep */ "v_pos", "comments", /*==*/ " " /* blank terminates */ }; + + const int efield_table_types[] = { 3, 2, 2, 2, 2, @@ -691,6 +693,18 @@ const char* const track_table_cols[] = " " /* blank terminates */ }; +const char* const dist_table_cols[] = +{ + "number", "x", "px", "y", "py", "t", "pt", + " " /* blank terminates */ +}; + +const int dist_table_types[] = +{ + 1, 2, 2, 2, 2, 2, 2 +}; + + const int track_table_cols_len = sizeof track_table_cols / sizeof track_table_cols[0]; const int tracksumm_table_types[] = diff --git a/src/mad_gcst.h b/src/mad_gcst.h index 5d0d9a59f..d07ef21d5 100644 --- a/src/mad_gcst.h +++ b/src/mad_gcst.h @@ -100,6 +100,9 @@ extern const char* const ptcnodetrack_table_cols[]; extern const int trackloss_table_types[]; extern const char* const trackloss_table_cols[]; +extern const char* const dist_table_cols[]; +extern const int dist_table_types[]; + extern const int dynap_table_types[]; extern const char* const dynap_table_cols[]; extern const int dynaptune_table_types[]; diff --git a/src/madx.h b/src/madx.h index 207239786..c76e3c20b 100644 --- a/src/madx.h +++ b/src/madx.h @@ -86,7 +86,8 @@ //#include "mad_elemrfc.h" // physics modules - +#include "../libs/DISTlib/source/distinterface.h" +#include "mad_dist.h" #include "mad_aper.h" #include "mad_dynap.h" #include "mad_emit.h" From 3194973162fd1537d1d0cb5659b0b9891e5e1815 Mon Sep 17 00:00:00 2001 From: Tobias Date: Fri, 20 Mar 2020 10:13:23 +0100 Subject: [PATCH 02/20] Not for the final implementation but to have it all stored. --- libs/DISTlib/CMakeLists.txt | 2 + libs/DISTlib/source/CMakeLists.txt | 6 + libs/DISTlib/source/distgeneration.c | 317 +++++ libs/DISTlib/source/distgeneration.h | 95 ++ libs/DISTlib/source/distinterface.c | 342 ++++++ libs/DISTlib/source/distinterface.h | 34 + libs/DISTlib/source/feigen.c | 1587 ++++++++++++++++++++++++++ libs/DISTlib/source/file_reader.c | 526 +++++++++ libs/DISTlib/source/file_reader.h | 18 + libs/DISTlib/source/helper.c | 148 +++ libs/DISTlib/source/helper.h | 17 + libs/DISTlib/source/outputdist.c | 33 + libs/DISTlib/source/outputdist.h | 10 + 13 files changed, 3135 insertions(+) create mode 100644 libs/DISTlib/CMakeLists.txt create mode 100644 libs/DISTlib/source/CMakeLists.txt create mode 100644 libs/DISTlib/source/distgeneration.c create mode 100644 libs/DISTlib/source/distgeneration.h create mode 100644 libs/DISTlib/source/distinterface.c create mode 100644 libs/DISTlib/source/distinterface.h create mode 100644 libs/DISTlib/source/feigen.c create mode 100644 libs/DISTlib/source/file_reader.c create mode 100644 libs/DISTlib/source/file_reader.h create mode 100644 libs/DISTlib/source/helper.c create mode 100644 libs/DISTlib/source/helper.h create mode 100644 libs/DISTlib/source/outputdist.c create mode 100644 libs/DISTlib/source/outputdist.h diff --git a/libs/DISTlib/CMakeLists.txt b/libs/DISTlib/CMakeLists.txt new file mode 100644 index 000000000..6cfd68090 --- /dev/null +++ b/libs/DISTlib/CMakeLists.txt @@ -0,0 +1,2 @@ +include(buildtypeCoverage) +add_subdirectory(source) diff --git a/libs/DISTlib/source/CMakeLists.txt b/libs/DISTlib/source/CMakeLists.txt new file mode 100644 index 000000000..ace06dcee --- /dev/null +++ b/libs/DISTlib/source/CMakeLists.txt @@ -0,0 +1,6 @@ +file(GLOB src_files *.c *.h) + +add_library(DISTlib STATIC ${src_files}) +if (NOT BUILD_SHARED_LIBS) + install(TARGETS DISTlib ARCHIVE DESTINATION lib) +endif () diff --git a/libs/DISTlib/source/distgeneration.c b/libs/DISTlib/source/distgeneration.c new file mode 100644 index 000000000..fcc9abf08 --- /dev/null +++ b/libs/DISTlib/source/distgeneration.c @@ -0,0 +1,317 @@ +#include +#include +#include +#include +#include "helper.h" +#include "distgeneration.h" +#include "distinterface.h" +#include "outputdist.h" +#include "file_reader.h" + + +void gensixcanonical(){ + + int counter = 0; + int type = dist->incoordtype; + double tc[dim]; + double normalized[dim], cancoord[dim]; + if(dist->ref->grid==1){ + generate_grid();//generate a grid + } + + for(int i =0; i< dist->totincoord; i++){ + for(int k=0; k<6; k++){ + tc[k]=dist->incoord[i]->coord[k]; + + } + + + if(type==0 || type==3){ + + action2normalized(tc, normalized); + normalized2canonical(normalized, cancoord); + } + else if(type==1){ + for(int k=0; k<6; k++){ + //normalized[k] = tc[k]; + normalized2canonical(tc, cancoord); + } + } + else if(type==2){ + for(int k=0; k<6; k++){ + cancoord[k]=tc[k]; + } + + } + + + if(particle_within_limits_physical(cancoord)==1 && particle_within_limits_normalized(normalized)){ + for(int p=0; poutcoord[counter]->coord[p] = cancoord[p]; + + } + dist->outcoord[i]->mass = dist->incoord[i]->mass; + dist->outcoord[i]->a = dist->incoord[i]->a; + dist->outcoord[i]->z = dist->incoord[i]->z; + counter++; + + } + } + dist->totoutcoord=counter; + dist->isDistrcalculated=1; + +} +int gettotalgridlength(){ + int totlength=1; + for(int i=0; iref->readinlength[i]; + } + dist->totincoord=totlength; + return totlength; +} +void generate_grid(){ + + int counter = 0; + double tc[6]; + double tmp[6]; + double tmp_n[6]; + double **readin; + int totallenght = gettotalgridlength(); + readin = (double**)malloc(dim*sizeof(double*)); + + for(int i=0; iref->readinlength[i]*sizeof(double)); + memcpy(readin[i],dist->incoord[i]->coord, dist->ref->readinlength[i]*sizeof(double)); + } + + deallocateincoord(); + + allocateincoord(totallenght); + + //dist->distout_normalized = (double**)malloc(getupperbound()*sizeof(double*)); + + for(int i =0; i< dist->ref->readinlength[0]; i++){ + for(int j =0; j< dist->ref->readinlength[1]; j++){ + for(int k =0; k< dist->ref->readinlength[2]; k++){ + for(int l =0; lref->readinlength[3]; l++){ + for(int m =0; m< dist->ref->readinlength[4]; m++){ + for(int n =0; n< dist->ref->readinlength[5]; n++){ + dist->incoord[counter]->coord[0]=readin[0][i]; + dist->incoord[counter]->coord[1]=readin[1][j]; + dist->incoord[counter]->coord[2]=readin[2][k]; + dist->incoord[counter]->coord[3]=readin[3][l]; + dist->incoord[counter]->coord[4]=readin[4][m]; + dist->incoord[counter]->coord[5]=readin[5][n]; + //dist->coord[5]->values[n]; + counter++; + + + } + } + } + } + } + } + dist->totincoord=counter; + dist->isDistrcalculated=1; +} + + +/*If emittance is defined it converts to canonical coordinates */ +void action2normalized(double acangl[6], double normalized[6]){ + + normalized[0]= sqrt(acangl[0])*cos(acangl[1]); + normalized[1]=-sqrt(acangl[0])*sin(acangl[1]); + normalized[2]= sqrt(acangl[2])*cos(acangl[3]); + normalized[3]=-sqrt(acangl[2])*sin(acangl[3]); + normalized[4]= sqrt(acangl[4])*cos(acangl[5]); + normalized[5]=-sqrt(acangl[4])*sin(acangl[5]); // used to devide with 1000 here before.. + +} + +void normalized2canonical(double normalized_in[6], double cancoord[6]){ + double normalized[6]; + normalized[0] = sqrt(dist->emitt->e1)*normalized_in[0]; + normalized[1] = sqrt(dist->emitt->e1)*normalized_in[1]; + normalized[2] = sqrt(dist->emitt->e2)*normalized_in[2]; + normalized[3] = sqrt(dist->emitt->e2)*normalized_in[3]; + normalized[4] = sqrt(dist->emitt->e3)*normalized_in[4]; + normalized[5] = sqrt(dist->emitt->e3)*normalized_in[5]; + + + if(dist->incoordtype==3) { + double lindp = 0; + double lindeltas=0; + double deltap = normalized[4]; + double deltas = normalized[5]; + double *xap; + double det = (dist->tas[4][4]*dist->tas[5][5] - dist->tas[4][5]*dist->tas[5][4]); + for(int i=0; i<4;i++){ + lindeltas = lindeltas+dist->tas[4][i]*normalized[i]; + lindp=lindp+dist->tas[5][i]*normalized[i]; + } + + lindp = deltap - lindp; + lindeltas = deltas - lindeltas; + + xap = (double*)malloc(2*sizeof(double)); + solve2by2eq(dist->tas[4][4], dist->tas[4][5], lindeltas, dist->tas[5][4], dist->tas[5][5], lindp, xap ); + normalized[4] = xap[0]; + normalized[5] = xap[1]; + + } + mtrx_vector_mult_pointer(dim,dim, dist->tas, normalized,cancoord); + +} + + +/*Checks if the particle is within the physical limit set by the user*/ +int particle_within_limits_physical(double *physical){ + + if(dist->cuts2apply->isset_p==0) return 1; + for(int i=0; icuts2apply->physical[i]->isset==1){ + if(physical[i] > dist->cuts2apply->physical[i]->min && physical[i] < dist->cuts2apply->physical[i]->max) return 0; + } + } + + return 1; + +} + +int particle_with_limits_action(int i, double value){ + + if(dist->cuts2apply->action[i]->isset==1){ + printf("iiiiiii %d \n", dist->cuts2apply->action[i]->isset); + if(value > pow(dist->cuts2apply->action[i]->min,2) && value < pow(dist->cuts2apply->action[i]->max,2)) return 1; + else return 0; + } + else{ + return 1; + } + + return 1; +} + +/*Checks if the particle is within the normalized limit set by the user*/ +int particle_within_limits_normalized(double *normalized){ + + if(dist->cuts2apply->isset_n==0) return 1; + for(int i=0; icuts2apply->normalized[i]->isset==1){ + if(normalized[i] > dist->cuts2apply->normalized[i]->min && normalized[i] < dist->cuts2apply->normalized[i]->max) return 0; + } + } + return 1; +} + +void createcoordinates(int index, double start, double stop, int length, int type){ + double temp [length]; + if(type ==0){ //Constant value + for(int i=0;i incoord[i]->coord[index] = start; + } + printf("tttttttttttttttt %d %f \n", index, start); + return; + } + + + else if(type==1){ //Linearly spaced intervalls + createLinearSpaced(length, start, stop, temp); + for(int i=0;i incoord[i]->coord[index] = temp[i]; + } + + } + else if(type==2){ //Exponentially spaced + createLinearSpaced(length, start, stop, temp); + for(int i=0;i incoord[i]->coord[index] = exp(temp[i]); + } + } + else if(type==3){ //Spaced with ^2 + createLinearSpaced(length, start, stop, temp); + for(int i=0;i incoord[i]->coord[index] = pow(temp[i],2); + + } + + } + else if(type==4){ // uniform random + for(int i=0;i incoord[i]->coord[index] = rand_uni(start, stop); + //printf("%f \n", dist->coord[index-1]->values[i] ); + } + } + + else if(type==5){ // Gaussian random (Here start is mean and stop is the standard deviation) + for(int i=0;i incoord[i]->coord[index] = randn(start, stop); + } + } + + else if(type==6){ // Rayleigh distribution + double tmp; + for(int i=0;i incoord[i]->coord[index] = tmp; + printf("%f \n", sqrt(tmp/2)); + } + } + else + issue_error("Unknown type of spacing"); +} + + +void createLinearSpaced(int length, double start, double stop, double *eqspaced ){ + + double distance = (stop-start)/length; + for(int i=0; i= 1 || W == 0); + + mult = sqrt ((-2 * log (W)) / W); + X1 = U1 * mult; + X2 = U2 * mult; + + call = !call; + + return (mu + sigma * (double) X1); +} + +double rand_uni(double low, double high) +{ + + return ( (double)rand() * ( high - low ) ) / (double)RAND_MAX + low; +} + +double randray(double mu, double sigma){ + double low = 0; + double high =1; + return pow((mu+sigma*sqrt((-2*log(rand_uni(low, high))))),2); +} \ No newline at end of file diff --git a/libs/DISTlib/source/distgeneration.h b/libs/DISTlib/source/distgeneration.h new file mode 100644 index 000000000..df6d6aa22 --- /dev/null +++ b/libs/DISTlib/source/distgeneration.h @@ -0,0 +1,95 @@ +struct distparam +{ + struct parameters** coord; + int disttype; + struct refparam* ref; + struct emittances* emitt; + double **tas; + double **invtas; + double *closedorbit; + int incoordtype; // This tells which type of coordinates the input is given. // 0-action angle, 1-Normalized, 2-physical, 3-mixed action angle and physical (lat) + struct coordinates** incoord; + struct coordinates** outcoord; + struct coordinates** gridin; + int totincoord; + int totoutcoord; + int isDistrcalculated; + int isallocated; + struct appliedcut* cuts2apply; +}; +struct refparam{ + int charge0; + int z0; + int a0; + double pc0; + double e0; + double mass0; + double beta0; + int en_like; + int time_like; + int ang_like; + struct emittances* emitt; + int *typeused; // This says which is used or a mixture of it.. (used for reading in and as a cross check) + // 0-action, 1-normalized, 2-physical + int *readinlength; + int grid; + +}; +struct coordinates +{ + double *coord; + double *readin; + double mass; + int charge; + int a; + int z; + +}; + +struct parameters +{ + double start; + double stop; + int length; + int type; //This gives the type of distribution, constant, linear, gaussian, + double * values; + int coordtype; +}; + +struct emittances{ + double e1, e2, e3; + double dp, deltas; +}; + +struct appliedcut{ + int isset_p; + int isset_n; + int isset_a; + struct cut** physical; + struct cut** normalized; + struct cut** action; +}; + +//void canonical2emittance_(double cancord[6], double emittance[3]); +//void dist2sixcoord_(); +//void action2canonical_(double acangl[6], double cancord[6], double acoord[6]); +//void action2sixinternal_(double tc[6], double *results, double *normalized); +//int checkdist(); +//void createrandomdist_(); +void gensixcanonical(); +int particle_within_limits_physical(double *physical); +int particle_within_limits_normalized(double *normalized); +int particle_with_limits_action(int i, double value); +void generatefromnormalized(); +void generatefrommixed(); +void generatefromaction(); +void generatefromphysical(); +void action2normalized(double acangl[6], double normalized[6]); +void normalized2canonical(double normalized[6], double cancoord[6]); +double randn(double mu, double sigma); +double rand_uni(double low, double high); +void createLinearSpaced(int length, double start, double stop, double *eqspaced ); +double randray(double mu, double sigma); +void createcoordinates(int index, double start, double stop, int length, int type); +int gettotalgridlength(); +void generate_grid(); \ No newline at end of file diff --git a/libs/DISTlib/source/distinterface.c b/libs/DISTlib/source/distinterface.c new file mode 100644 index 000000000..bb277618a --- /dev/null +++ b/libs/DISTlib/source/distinterface.c @@ -0,0 +1,342 @@ +#include +#include +#include +#include +#include "helper.h" +#include "distinterface.h" +#include "distgeneration.h" +#include "outputdist.h" +#include "file_reader.h" + + +/* +This allocates the the memory for the distributions +*/ +void initializedistribution(int numberOfDist){ + dist = (struct distparam*)malloc((numberOfDist)*sizeof(struct distparam)); + dim = 6; + + for(int i = 0; i ref = (struct refparam*)malloc(sizeof(struct refparam)); + (dist + i)->coord = (struct parameters**)malloc(dim*sizeof(struct parameters*)); + (dist + i)->emitt = (struct emittances*)malloc(sizeof(struct emittances)); + (dist + i)->cuts2apply = (struct appliedcut*)malloc(sizeof(struct appliedcut)); + (dist + i)->cuts2apply->physical = (struct cut**)malloc(dim*sizeof(struct cut*)); + (dist + i)->cuts2apply->normalized = (struct cut**)malloc(dim*sizeof(struct cut*)); + (dist + i)->cuts2apply->action = (struct cut**)malloc(dim*sizeof(struct cut*)); + (dist + i)->tas = (double**)malloc(dim*sizeof(double*)); + (dist + i)->invtas = (double**)malloc(dim*sizeof(double*)); + (dist + i)->closedorbit = (double*)malloc(dim*sizeof(double)); + (dist + i)->isDistrcalculated = 0; + (dist + i)->ref->e0=0; + (dist + i)->ref->pc0=0; + (dist + i)->ref->a0=1; + (dist + i)->ref->z0=1; + (dist + i)->ref->mass0=0; + (dist + i)->ref->charge0=1; + (dist + i)->ref->en_like=-1; + (dist + i)->ref->time_like=-1; + (dist + i)->ref->ang_like=-1; + (dist + i)->totincoord=-1; + (dist + i)->ref->typeused = (int*)malloc(dim*sizeof(int)); + (dist + i)->ref->readinlength = (int*)malloc(dim*sizeof(int)); + (dist + i)->ref->grid=-1; + for(int k=0; ktas[k] =(double*)malloc(dim*sizeof(double)); + (dist + i)->invtas[k] =(double*)malloc(dim*sizeof(double)); + } + (dist + i)->incoordtype =-1; + (dist + i)->disttype = 0; + for(int j=0; jcuts2apply->physical[j] = (struct cut*)malloc(sizeof(struct cut)); + (dist + i)->cuts2apply->normalized[j] = (struct cut*)malloc(sizeof(struct cut)); + (dist + i)->cuts2apply->action[j] = (struct cut*)malloc(sizeof(struct cut)); + (dist +i)->coord[j] = (struct parameters*)malloc(sizeof(struct parameters)); + (dist +i)->coord[j]->start=0; + (dist +i)->coord[j]->stop=0; + (dist +i)->coord[j]->length=1; + (dist +i)->coord[j]->type=0; + (dist +i)->closedorbit[j]=0; + } + for (int i = 0; i < 6; i++) + { + for (int k = 0; k < 6; k++) + { + dist->tas [i][k] = 0; + dist->invtas[i][k] = 0; + + } + } + } + diststart=dist; + +} + +void sete0andmass0(double energy0, double mass0){ + dist->ref->mass0 = mass0; + dist->ref->e0 = energy0; + dist->ref->pc0 = energy2momentum(dist->ref->e0,dist->ref->mass0); + dist->ref->beta0 = (dist->ref->pc0)/(dist->ref->e0); +} + +void setdistribution(int ndist){ + dist = diststart + ndist; +} + +void setemitt12(double e1, double e2){ + dist->emitt->e1=e1; + dist->emitt->e2=e2; +} + +void setemitt3(double e3){ + dist->emitt->e3=e3; +} + +void settasmatrix(double *tas){ + for(int i =0; itas[i][j] = tas[j+i*dim]; + } + } +} + +void settasmatrixtranspose(double *tas){ + for(int i =0; itas[j][i] = tas[j+i*dim]; + } + } +} + +void setactionanglecut(int variable, double min, double max){ + dist->cuts2apply->isset_a=1; + dist->cuts2apply->action[variable]->min=min; + dist->cuts2apply->action[variable]->max=max; + dist->cuts2apply->action[variable]->isset=1; +} + +void settasmatrix_element(double value, int row, int column){ + + dist->tas[row][column] = value; + +} +/*Create a TAS matrix using E-T, assuming uncoupled system */ +void settwisstas(double betax, double alfax, double betay, double alfay){ + + dist->tas[0][0] = sqrt(betax); + dist->tas[1][0] =-(alfax)/sqrt(betax); + dist->tas[1][1] =-1/sqrt(betax); + + + dist->tas[2][2] = sqrt(betay); + dist->tas[3][2] =-alfay/sqrt(betay); + dist->tas[3][3] =-1/sqrt(betay); + + dist->tas[5][5] =1; // This is not a matched longitudinal but enables to change the energy + dist->tas[4][4] =1; + +} + +void setdisptas(double dx, double dpx, double dy, double dpy){ + + dist->tas[0][5] = dx; + dist->tas[1][5] = dpx; + dist->tas[2][5] = dy; + dist->tas[3][5] = dpy; + +} + + +// 0 -action angle, 1- normalized coordinates, 2-physical, 3-mixed +void setcoords(double *xn, double *xpn, double *yn, double *ypn, double *zn, double *zpn, int totparticles, int coordtype){ + allocateincoord(totparticles); + for(int i=0; iincoord[i]->coord[0] = xn[i]; + dist->incoord[i]->coord[1] = xpn[i]; + dist->incoord[i]->coord[2] = yn[i]; + dist->incoord[i]->coord[3] = ypn[i]; + dist->incoord[i]->coord[4] = zn[i]; + dist->incoord[i]->coord[5] = zpn[i]; + } + + dist->totincoord = totparticles; + dist->incoordtype = coordtype; +} + +void settotalsteps(int totgenerate){ + dist->totincoord = totgenerate; +} + +void setscan_para_diagonal(int variable, int variable_type, int type, double start, double stop){ + if(dist->totincoord==-1) + issue_error("Must set a total steps before you can set a parameter"); + + if(dist->isallocated!=1){ + allocateincoord(dist->totincoord); + } + dist->ref->typeused[variable] = variable_type; + dist->incoordtype=variable_type; // this is wrong + + createcoordinates(variable, start, stop, dist->totincoord,type); +} + +void setscan_para_grid(int variable,int variable_type,int type, double start, double stop, int length){ + dist->ref->readinlength[variable] = length; + + if(length>MAX_LENGTH) + issue_error("For grid scans you have requested to many particles."); + if(dist->isallocated!=1) + allocateincoord(MAX_LENGTH); + + dist->ref->typeused[variable] = variable_type; + dist->incoordtype=variable_type; + createcoordinates(variable, start, stop, dist->totincoord,type); + dist->ref->grid=1; + +} + +void addclosedorbit(double *clo){ + for(int i=0; iclosedorbit[i] = clo[i]; + } +} + + +void setphysicalcut(int variable, double min, double max){ + dist->cuts2apply->isset_p=1; + dist->cuts2apply->physical[variable-1]->min=min; + dist->cuts2apply->physical[variable-1]->max=max; + dist->cuts2apply->physical[variable-1]->isset=1; + +} + +void setnormalizedcut(int variable, double min, double max){ + dist->cuts2apply->isset_n=1; + dist->cuts2apply->normalized[variable-1]->min=min; + dist->cuts2apply->normalized[variable-1]->max=max; + dist->cuts2apply->normalized[variable-1]->isset=1; + +} +void getarraylength(int *totlength){ + if(dist->isDistrcalculated ==0){ + + gensixcanonical(); + } + *totlength=dist->totoutcoord; +} + +void get6trackcoord(double *x, double *xp, double *y, double *yp, double *sigma, double *deltap, int *totparticles){ + double tmp[6]; + int nparticles; + if(dist->isDistrcalculated ==0){ + gensixcanonical(); + + } + if(dist->totoutcoord < *totparticles) + nparticles = dist->totoutcoord; + else + nparticles = *totparticles; + + for(int i=0; i < nparticles; i++){ + canonical2six(dist->outcoord[i]->coord, dist->ref->beta0, dist->ref->pc0, dist->ref->mass0, dist->incoord[i]->mass, tmp); + x[i] = tmp[0]; + xp[i] = tmp[1]; + y[i] = tmp[2]; + yp[i] = tmp[3]; + sigma[i] = tmp[4]; + deltap[i] = tmp[5]; + } + *totparticles=nparticles; +} + +void getunconvertedcoord(double *x, double *xp, double *y, double *yp, double *sigma, double *deltap, int *totparticles){ + double tmp[6]; + int nparticles; + if(dist->isDistrcalculated ==0){ + gensixcanonical(); + + } + if(dist->totoutcoord < *totparticles) + nparticles = dist->totoutcoord; + else + nparticles = *totparticles; + + for(int i=0; i < nparticles; i++){ + x[i] = dist->outcoord[i]->coord[0]; + xp[i] = dist->outcoord[i]->coord[1]; + y[i] = dist->outcoord[i]->coord[2]; + yp[i] = dist->outcoord[i]->coord[3]; + sigma[i] = dist->outcoord[i]->coord[4]; + deltap[i] = dist->outcoord[i]->coord[5]; + } + *totparticles=nparticles; +} + +void readtasmatrixfile(const char* filename_in){ + FILE * fp; + char line [500]; /* or other suitable maximum line size */ + + fp = fopen(filename_in, "r+"); + int i=0; + if(fp==NULL){ + printf("No such file\n"); + exit(1); + } + while ( fgets ( line, sizeof line, fp ) != NULL ){ /* read a line */ + printf("%s \n", line); + sscanf(line, "%lf %lf %lf %lf %lf %lf", &dist->tas[i][1], &dist->tas[i][1], &dist->tas[i][2], &dist->tas[i][3], &dist->tas[i][4], &dist->tas[i][5]); + //printf("%f %f %f %f %f %f \n", dist->tas[i][1], dist->tas[i][1], dist->tas[i][2], dist->tas[i][3], dist->tas[i][4], dist->tas[i][5]); + i++; + } + fclose(fp); +} + + +void getrefpara(double *energy0, double *mass0, int *a0, int *z0){ + *energy0=dist->ref->e0; + *mass0=dist->ref->mass0; + *a0=dist->ref->a0; + *z0=dist->ref->z0; +} +int readfile_f(const char* filename_in, int strlen){ + char filename [strlen]; + strncpy(filename, filename_in, strlen); + filename[strlen] = '\0'; + readfile(filename); +} + +int writefile_f(const char* filename_in, int strlen){ + char filename [strlen]; + strncpy(filename, filename_in, strlen); + filename[strlen] = '\0'; + print2file(filename); +} + +void free_distribution(int i){ + + free((dist+i)->closedorbit); + free((dist+i)->tas); + free((dist+i)->invtas); + + for(int j=0; jcuts2apply->physical[j]) ; + free((dist + i)->cuts2apply->normalized[j]) ; + free((dist + i)->cuts2apply->action[j]) ; + free((dist +i)->coord[j]); + } + free((dist+i)->coord); + free((dist + i)->cuts2apply->physical); + free((dist + i)->cuts2apply->normalized); + free((dist + i)->cuts2apply->action); + free((dist + i)->cuts2apply); + + free((dist + i)->ref); + free((dist + i)->emitt); + free(dist+i); + + +} \ No newline at end of file diff --git a/libs/DISTlib/source/distinterface.h b/libs/DISTlib/source/distinterface.h new file mode 100644 index 000000000..d2a9756ca --- /dev/null +++ b/libs/DISTlib/source/distinterface.h @@ -0,0 +1,34 @@ +struct cut{ + int isset; + double min; + double max; +}; + + +struct distparam* dist; +struct distparam* diststart; +int dim; +int distn; +int writefile_f(const char* filename_in, int strlen); +int readfile_f(const char* filename_in, int strlen); +void readtasmatrixfile(const char* filename_in); +void initializedistribution(int numberOfDist); +void setdistribution(int ndist); +void setphysicalcut(int variable, double min, double max); +void setnormalizedcut(int variable, double min, double max); + +void settotalsteps(int totgenerate); +void sete0andmass0(double energy0, double mass0); +void setemitt12(double e1, double e2); +void setemitt3(double e3); +void settasmatrix(double * tas); +void settasmatrixtranspose(double *tas); +void setcoords(double *xn, double *xnp, double *yn, double *ynp, double *zn, double *znp, int totparticles, int coordtype); +void get6trackcoord(double *x, double *xp, double *y, double *yp, double *sigma, double *deltap, int *totparticles); +void setscan_para_diagonal(int variable, int variable_type, int type, double start, double stop); +void setscan_para_grid(int variable,int variable_type,int type, double start, double stop, int length); +void setdisptas(double dx, double dpx, double dy, double dpy); +void settwisstas(double betax, double alfax, double betay, double alfay); +void getunconvertedcoord(double *x, double *xp, double *y, double *yp, double *sigma, double *deltap, int *totparticles); +void setactionanglecut(int variable, double min, double max); +void free_distribution(int i); \ No newline at end of file diff --git a/libs/DISTlib/source/feigen.c b/libs/DISTlib/source/feigen.c new file mode 100644 index 000000000..a7eb14b85 --- /dev/null +++ b/libs/DISTlib/source/feigen.c @@ -0,0 +1,1587 @@ +/* + This source file is adapted from feigen.c that comes with the book + Numeric Algorithm with C by Frank Uhlig et al. I cannot find the + license of the original source codes. I release my modifications under + the MIT license. The modified version requires C99 as it uses complex + numbers. I may modify the code if this is a concern. + */ + +/* The MIT License + + Copyright (c) 1996 Frank Uhlig et al. + 2009 Genome Research Ltd (GRL). + + Permission is hereby granted, free of charge, to any person obtaining + a copy of this software and associated documentation files (the + "Software"), to deal in the Software without restriction, including + without limitation the rights to use, copy, modify, merge, publish, + distribute, sublicense, and/or sell copies of the Software, and to + permit persons to whom the Software is furnished to do so, subject to + the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS + BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN + ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ + +/* Contact: Heng Li */ + + +/*.FE{C 7.8} + {QR Algorithm} + {Eigenvalues and Eigenvectors of a Matrix via the + QR Algorithm}*/ + +/*.BE*/ +/* ------------------------- MODULE feigen.c ------------------------ */ + +#include +#include /* for DBL_EPSILON */ +#include +#include +#include + +#define REAL double +#define TRUE 1 +#define FALSE 0 +#define ZERO 0. +#define ONE 1. +#define TWO 2. +#define VEKTOR 0 +#define MACH_EPS DBL_EPSILON + +#define ABS(x) (((x) >= 0.)? (x) : -(x)) +#define SQRT(x) sqrt(x) +#define SQR(x) ((x) * (x)) +#define SWAP(typ, a, b) { typ t; t = (a); (a) = (b); (b) = t; } +#define BASIS basis() + +typedef struct { + int n, max; + REAL *mem; +} vmblock_t; + +static void *vminit(void) +{ + return (vmblock_t*)calloc(1, sizeof(vmblock_t)); +} +static int vmcomplete(void *vmblock) +{ + return 1; +} +static void vmfree(void *vmblock) +{ + vmblock_t *vmb = (vmblock_t*)vmblock; + free(vmb->mem); free(vmblock); +} +static void *vmalloc(void *vmblock, int typ, size_t zeilen, size_t spalten) +{ + vmblock_t *vmb = (vmblock_t*)vmblock; + double *ret = 0; + if (typ == 0) { + if (vmb->n + zeilen > vmb->max) { + vmb->max = vmb->n + zeilen; + vmb->mem = (REAL*)realloc(vmb->mem, vmb->max * sizeof(REAL)); + } + ret = vmb->mem + vmb->n; + vmb->n += zeilen; + } + return ret; +} + +static int basis(void) /* find basis used for computer number representation */ +/*.IX{basis}*/ + +/*********************************************************************** +* Find the basis for representing numbers on the computer, if not * +* already done. * +* * +* global names used: * +* ================== * +* REAL, ONE, TWO * +***********************************************************************/ + +{ + REAL x, + eins, + b; + + x = eins = b = ONE; + + while ((x + eins) - x == eins) + x *= TWO; + while ((x + b) == x) + b *= TWO; + + + return (int)((x + b) - x); +} + + +#define MAXIT 50 /* Maximal number of */ + /* iterations per eigenvalue */ + + +/*--------------------------------------------------------------------* + * Aux functions for eigen * + *--------------------------------------------------------------------*/ + + +static int balance /* balance a matrix .........................*/ +/*.IX{balance}*/ + (int n, /* size of matrix ..............*/ + REAL * mat[], /* matrix ......................*/ + REAL scal[], /* Scaling data ................*/ + int * low, /* first relevant row index ....*/ + int * high, /* last relevant row index .....*/ + int basis /* base of computer numbers ....*/ + ) +/*====================================================================* + * * + * balance balances the matrix mat so that the rows with zero entries* + * off the diagonal are isolated and the remaining columns and rows * + * are resized to have one norm close to 1. * + * * + *====================================================================* + * * + * Input parameters: * + * ================ * + * n int n; ( n > 0 ) * + * Dimension of mat * + * mat REAL *mat[n]; * + * n x n input matrix * + * basis int basis; * + * Base of number representaion in the given computer * + * (see BASIS) * + * * + * Output parameters: * + * ================== * + * mat REAL *mat[n]; * + * scaled matrix * + * low int *low; * + * high int *high; * + * the rows 0 to low-1 and those from high to n-1 * + * contain isolated eigenvalues (only nonzero entry on * + * the diagonal) * + * scal REAL scal[]; * + * the vector scal contains the isolated eigenvalues in * + * the positions 0 to low-1 and high to n-1, its other * + * components contain the scaling factors for * + * transforming mat. * + * * + *====================================================================* + * * + * Macros: SWAP, ABS * + * ======= * + * * + *====================================================================* + * * + * Constants used : TRUE, FALSE * + * =================== * + * * + *====================================================================*/ +{ + register int i, j; + int iter, k, m; + REAL b2, r, c, f, g, s; + + b2 = (REAL) (basis * basis); + m = 0; + k = n - 1; + + do + { + iter = FALSE; + for (j = k; j >= 0; j--) + { + for (r = ZERO, i = 0; i <= k; i++) + if (i != j) r += ABS (mat[j][i]); + + if (r == ZERO) + { + scal[k] = (REAL) j; + if (j != k) + { + for (i = 0; i <= k; i++) SWAP (REAL, mat[i][j], mat[i][k]) + for (i = m; i < n; i++) SWAP (REAL, mat[j][i], mat[k][i]) + } + k--; + iter = TRUE; + } + } /* end of j */ + } /* end of do */ + while (iter); + + do + { + iter = FALSE; + for (j = m; j <= k; j++) + { + for (c = ZERO, i = m; i <= k; i++) + if (i != j) c += ABS (mat[i][j]); + if (c == ZERO) + { + scal[m] = (REAL) j; + if ( j != m ) + { + for (i = 0; i <= k; i++) SWAP (REAL, mat[i][j], mat[i][m]) + for (i = m; i < n; i++) SWAP (REAL, mat[j][i], mat[m][i]) + } + m++; + iter = TRUE; + } + } /* end of j */ + } /* end of do */ + while (iter); + + *low = m; + *high = k; + for (i = m; i <= k; i++) scal[i] = ONE; + + do + { + iter = FALSE; + for (i = m; i <= k; i++) + { + for (c = r = ZERO, j = m; j <= k; j++) + if (j !=i) + { + c += ABS (mat[j][i]); + r += ABS (mat[i][j]); + } + g = r / basis; + f = ONE; + s = c + r; + + while (c < g) + { + f *= basis; + c *= b2; + } + + g = r * basis; + while (c >= g) + { + f /= basis; + c /= b2; + } + + if ((c + r) / f < (REAL)0.95 * s) + { + g = ONE / f; + scal[i] *= f; + iter = TRUE; + for (j = m; j < n; j++ ) mat[i][j] *= g; + for (j = 0; j <= k; j++ ) mat[j][i] *= f; + } + } + } + while (iter); + + return (0); +} + + +static int balback /* reverse balancing ........................*/ +/*.IX{balback}*/ + (int n, /* Dimension of matrix .........*/ + int low, /* first nonzero row ...........*/ + int high, /* last nonzero row ............*/ + REAL scal[], /* Scaling data ................*/ + REAL * eivec[] /* Eigenvectors ................*/ + ) +/*====================================================================* + * * + * balback reverses the balancing of balance for the eigenvactors. * + * * + *====================================================================* + * * + * Input parameters: * + * ================ * + * n int n; ( n > 0 ) * + * Dimension of mat * + * low int low; * + * high int high; see balance * + * eivec REAL *eivec[n]; * + * Matrix of eigenvectors, as computed in qr2 * + * scal REAL scal[]; * + * Scaling data from balance * + * * + * Output parameter: * + * ================ * + * eivec REAL *eivec[n]; * + * Non-normalized eigenvectors of the original matrix * + * * + * Macros : SWAP() * + * ======== * + * * + *====================================================================*/ +{ + register int i, j, k; + REAL s; + + for (i = low; i <= high; i++) + { + s = scal[i]; + for (j = 0; j < n; j++) eivec[i][j] *= s; + } + + for (i = low - 1; i >= 0; i--) + { + k = (int) scal[i]; + if (k != i) + for (j = 0; j < n; j++) SWAP (REAL, eivec[i][j], eivec[k][j]) + } + + for (i = high + 1; i < n; i++) + { + k = (int) scal[i]; + if (k != i) + for (j = 0; j < n; j++) SWAP (REAL, eivec[i][j], eivec[k][j]) + } + return (0); +} + + + +static int elmhes /* reduce matrix to upper Hessenberg form ....*/ +/*.IX{elmhes}*/ + (int n, /* Dimension of matrix .........*/ + int low, /* first nonzero row ...........*/ + int high, /* last nonzero row ............*/ + REAL * mat[], /* input/output matrix .........*/ + int perm[] /* Permutation vector ..........*/ + ) +/*====================================================================* + * * + * elmhes transforms the matrix mat to upper Hessenberg form. * + * * + *====================================================================* + * * + * Input parameters: * + * ================ * + * n int n; ( n > 0 ) * + * Dimension of mat * + * low int low; * + * high int high; see balance * + * mat REAL *mat[n]; * + * n x n matrix * + * * + * Output parameter: * + * ================= * + * mat REAL *mat[n]; * + * upper Hessenberg matrix; additional information on * + * the transformation is stored in the lower triangle * + * perm int perm[]; * + * Permutation vector for elmtrans * + * * + *====================================================================* + * * + * Macros: SWAP, ABS * + * ======= * + * * + *====================================================================*/ +{ + register int i, j, m; + REAL x, y; + + for (m = low + 1; m < high; m++) + { + i = m; + x = ZERO; + for (j = m; j <= high; j++) + if (ABS (mat[j][m-1]) > ABS (x)) + { + x = mat[j][m-1]; + i = j; + } + + perm[m] = i; + if (i != m) + { + for (j = m - 1; j < n; j++) SWAP (REAL, mat[i][j], mat[m][j]) + for (j = 0; j <= high; j++) SWAP (REAL, mat[j][i], mat[j][m]) + } + + if (x != ZERO) + { + for (i = m + 1; i <= high; i++) + { + y = mat[i][m-1]; + if (y != ZERO) + { + y = mat[i][m-1] = y / x; + for (j = m; j < n; j++) mat[i][j] -= y * mat[m][j]; + for (j = 0; j <= high; j++) mat[j][m] += y * mat[j][i]; + } + } /* end i */ + } + } /* end m */ + + return (0); +} + + +static int elmtrans /* copy to Hessenberg form .................*/ +/*.IX{elmtrans}*/ + (int n, /* Dimension of matrix .........*/ + int low, /* first nonzero row ...........*/ + int high, /* last nonzero row ............*/ + REAL * mat[], /* input matrix ................*/ + int perm[], /* row permutations ............*/ + REAL * h[] /* Hessenberg matrix ...........*/ + ) +/*====================================================================* + * * + * elmtrans copies the Hessenberg matrix stored in mat to h. * + * * + *====================================================================* + * * + * Input parameters: * + * ================ * + * n int n; ( n > 0 ) * + * Dimension of mat and eivec * + * low int low; * + * high int high; see balance * + * mat REAL *mat[n]; * + * n x n input matrix * + * perm int *perm; * + * Permutation data from elmhes * + * * + * Output parameter: * + * ================ * + * h REAL *h[n]; * + * Hessenberg matrix * + * * + *====================================================================*/ +{ + register int k, i, j; + + for (i = 0; i < n; i++) + { + for (k = 0; k < n; k++) h[i][k] = ZERO; + h[i][i] = ONE; + } + + for (i = high - 1; i > low; i--) + { + j = perm[i]; + for (k = i + 1; k <= high; k++) h[k][i] = mat[k][i-1]; + if ( i != j ) + { + for (k = i; k <= high; k++) + { + h[i][k] = h[j][k]; + h[j][k] = ZERO; + } + h[j][i] = ONE; + } + } + + return (0); +} + + +/* ------------------------------------------------------------------ */ + +static int orthes /* reduce orthogonally to upper Hessenberg form */ +/*.IX{orthes}*/ + ( + int n, /* Dimension of matrix */ + int low, /* [low,low]..[high,high]: */ + int high, /* submatrix to be reduced */ + REAL *mat[], /* input/output matrix */ + REAL d[] /* reduction information */ + ) /* error code */ + +/*********************************************************************** +* This function reduces matrix mat to upper Hessenberg form by * +* Householder transformations. All details of the transformations are * +* stored in the remaining triangle of the Hessenberg matrix and in * +* vector d. * +* * +* Input parameters: * +* ================= * +* n dimension of mat * +* low \ rows 0 to low-1 and high+1 to n-1 contain isolated * +* high > eigenvalues, i. e. eigenvalues corresponding to * +* / eigenvectors that are multiples of unit vectors * +* mat [0..n-1,0..n-1] matrix to be reduced * +* * +* Output parameters: * +* ================== * +* mat the desired Hessenberg matrix together with the first part * +* of the reduction information below the subdiagonal * +* d [0..n-1] vector with the remaining reduction information * +* * +* Return value: * +* ============= * +* Error code. This can only be the value 0 here. * +* * +* global names used: * +* ================== * +* REAL, MACH_EPS, ZERO, SQRT * +* * +************************************************************************ +* Literature: Numerical Mathematics 12 (1968), pages 359 and 360 * +***********************************************************************/ + +{ + int i, j, m; /* loop variables */ + REAL s, /* Euclidian norm sigma of the subdiagonal column */ + /* vector v of mat, that shall be reflected into a */ + /* multiple of the unit vector e1 = (1,0,...,0) */ + /* (v = (v1,..,v(high-m+1)) */ + x = ZERO, /* first element of v in the beginning, then */ + /* summation variable in the actual Householder */ + /* transformation */ + y, /* sigma^2 in the beginning, then ||u||^2, with */ + /* u := v +- sigma * e1 */ + eps; /* tolerance for checking if the transformation is */ + /* valid */ + + eps = (REAL)128.0 * MACH_EPS; + + for (m = low + 1; m < high; m++) + { + for (y = ZERO, i = high; i >= m; i--) + x = mat[i][m - 1], + d[i] = x, + y = y + x * x; + if (y <= eps) + s = ZERO; + else + { + s = (x >= ZERO) ? -SQRT(y) : SQRT(y); + y -= x * s; + d[m] = x - s; + + for (j = m; j < n; j++) /* multiply mat from the */ + { /* left by (E-(u*uT)/y) */ + for (x = ZERO, i = high; i >= m; i--) + x += d[i] * mat[i][j]; + for (x /= y, i = m; i <= high; i++) + mat[i][j] -= x * d[i]; + } + + for (i = 0; i <= high; i++) /* multiply mat from the */ + { /* right by (E-(u*uT)/y) */ + for (x = ZERO, j = high; j >= m; j--) + x += d[j] * mat[i][j]; + for (x /= y, j = m; j <= high; j++) + mat[i][j] -= x * d[j]; + } + } + + mat[m][m - 1] = s; + } + + return 0; + +} /* --------------------------- orthes -------------------------- */ + + +/* ------------------------------------------------------------------ */ + +static int orttrans /* compute orthogonal transformation matrix */ +/*.IX{orttrans}*/ + ( + int n, /* Dimension of matrix */ + int low, /* [low,low]..[high,high]: submatrix */ + int high, /* affected by the reduction */ + REAL *mat[], /* Hessenberg matrix, reduction inf. */ + REAL d[], /* remaining reduction information */ + REAL *v[] /* transformation matrix */ + ) /* error code */ + +/*********************************************************************** +* compute the matrix v of accumulated transformations from the * +* information left by the Householder reduction of matrix mat to upper * +* Hessenberg form below the Hessenberg matrix in mat and in the * +* vector d. The contents of the latter are destroyed. * +* * +* Input parameters: * +* ================= * +* n dimension of mat * +* low \ rows 0 to low-1 and high+1 to n-1 contain isolated * +* high > eigenvalues, i. e. eigenvalues corresponding to * +* / eigenvectors that are multiples of unit vectors * +* mat [0..n-1,0..n-1] matrix produced by `orthes' giving the * +* upper Hessenberg matrix and part of the information on the * +* orthogonal reduction * +* d [0..n-1] vector with the remaining information on the * +* orthogonal reduction to upper Hessenberg form * +* * +* Output parameters: * +* ================== * +* d input vector destroyed by this function * +* v [0..n-1,0..n-1] matrix defining the similarity reduction * +* to upper Hessenberg form * +* * +* Return value: * +* ============= * +* Error code. This can only be the value 0 here. * +* * +* benutzte globale Namen: * +* ======================= * +* REAL, ZERO, ONE * +* * +************************************************************************ +* Literature: Numerical Mathematics 16 (1970), page 191 * +***********************************************************************/ + +{ + int i, j, m; /* loop variables */ + REAL x, /* summation variable in the */ + /* Householder transformation */ + y; /* sigma respectively */ + /* sigma * (v1 +- sigma) */ + + for (i = 0; i < n; i++) /* form the unit matrix in v */ + { + for (j = 0; j < n; j++) + v[i][j] = ZERO; + v[i][i] = ONE; + } + + for (m = high - 1; m > low; m--) /* apply the transformations */ + { /* that reduced mat to upper */ + y = mat[m][m - 1]; /* Hessenberg form also to the */ + /* unit matrix in v. This */ + if (y != ZERO) /* produces the desired */ + { /* transformation matrix in v. */ + y *= d[m]; + for (i = m + 1; i <= high; i++) + d[i] = mat[i][m - 1]; + for (j = m; j <= high; j++) + { + for (x = ZERO, i = m; i <= high; i++) + x += d[i] * v[i][j]; + for (x /= y, i = m; i <= high; i++) + v[i][j] += x * d[i]; + } + } + } + + return 0; + +} /* -------------------------- orttrans ------------------------- */ + + +static int hqrvec /* compute eigenvectors ......................*/ +/*.IX{hqrvec}*/ + (int n, /* Dimension of matrix .......*/ + int low, /* first nonzero row .........*/ + int high, /* last nonzero row ..........*/ + REAL * h[], /* upper Hessenberg matrix ...*/ + REAL wr[], /* Real parts of evalues .....*/ + REAL wi[], /* Imaginary parts of evalues */ + REAL * eivec[] /* Eigenvectors ..............*/ + ) +/*====================================================================* + * * + * hqrvec computes the eigenvectors for the eigenvalues found in hqr2* + * * + *====================================================================* + * * + * Input parameters: * + * ================ * + * n int n; ( n > 0 ) * + * Dimension of mat and eivec, number of eigenvalues. * + * low int low; * + * high int high; see balance * + * h REAL *h[n]; * + * upper Hessenberg matrix * + * wr REAL wr[n]; * + * Real parts of the n eigenvalues. * + * wi REAL wi[n]; * + * Imaginary parts of the n eigenvalues. * + * * + * Output parameter: * + * ================ * + * eivec REAL *eivec[n]; * + * Matrix, whose columns are the eigenvectors * + * * + * Return value : * + * ============= * + * = 0 all ok * + * = 1 h is the zero matrix. * + * * + *====================================================================* + * * + * function in use : * + * ================== * + * * + * int comdiv(): complex division * + * * + *====================================================================* + * * + * Constants used : MACH_EPS * + * ================= * + * * + * Macros : SQR, ABS * + * ======== * + * * + *====================================================================*/ +{ + int i, j, k; + int l, m, en, na; + REAL p, q, r = ZERO, s = ZERO, t, w, x, y, z = ZERO, + ra, sa, vr, vi, norm; + + for (norm = ZERO, i = 0; i < n; i++) /* find norm of h */ + { + for (j = i; j < n; j++) norm += ABS(h[i][j]); + } + + if (norm == ZERO) return (1); /* zero matrix */ + + for (en = n - 1; en >= 0; en--) /* transform back */ + { + p = wr[en]; + q = wi[en]; + na = en - 1; + if (q == ZERO) + { + m = en; + h[en][en] = ONE; + for (i = na; i >= 0; i--) + { + w = h[i][i] - p; + r = h[i][en]; + for (j = m; j <= na; j++) r += h[i][j] * h[j][en]; + if (wi[i] < ZERO) + { + z = w; + s = r; + } + else + { + m = i; + if (wi[i] == ZERO) + h[i][en] = -r / ((w != ZERO) ? (w) : (MACH_EPS * norm)); + else + { /* Solve the linear system: */ + /* | w x | | h[i][en] | | -r | */ + /* | | | | = | | */ + /* | y z | | h[i+1][en] | | -s | */ + + x = h[i][i+1]; + y = h[i+1][i]; + q = SQR (wr[i] - p) + SQR (wi[i]); + h[i][en] = t = (x * s - z * r) / q; + h[i+1][en] = ( (ABS(x) > ABS(z) ) ? + (-r -w * t) / x : (-s -y * t) / z); + } + } /* wi[i] >= 0 */ + } /* end i */ + } /* end q = 0 */ + + else if (q < ZERO) + { + m = na; + if (ABS(h[en][na]) > ABS(h[na][en])) + { + h[na][na] = - (h[en][en] - p) / h[en][na]; + h[na][en] = - q / h[en][na]; + } else { + /* comdiv(-h[na][en], ZERO, h[na][na]-p, q, &h[na][na], &h[na][en]); */ + REAL complex c; + c = -h[na][en] / (h[na][na]-p + q * I); + h[na][na] = creal(c); h[na][en] = cimag(c); + } + + h[en][na] = ONE; + h[en][en] = ZERO; + for (i = na - 1; i >= 0; i--) + { + w = h[i][i] - p; + ra = h[i][en]; + sa = ZERO; + for (j = m; j <= na; j++) + { + ra += h[i][j] * h[j][na]; + sa += h[i][j] * h[j][en]; + } + + if (wi[i] < ZERO) + { + z = w; + r = ra; + s = sa; + } + else + { + m = i; + if (wi[i] == ZERO) { + /* comdiv (-ra, -sa, w, q, &h[i][na], &h[i][en]); */ + REAL complex c; + c = (-ra - sa * I) / (w + q * I); + h[i][na] = creal(c); h[i][en] = cimag(c); + } else + { + + /* solve complex linear system: */ + /* | w+i*q x | | h[i][na] + i*h[i][en] | | -ra+i*sa | */ + /* | | | | = | | */ + /* | y z+i*q| | h[i+1][na]+i*h[i+1][en]| | -r+i*s | */ + + x = h[i][i+1]; + y = h[i+1][i]; + vr = SQR (wr[i] - p) + SQR (wi[i]) - SQR (q); + vi = TWO * q * (wr[i] - p); + if (vr == ZERO && vi == ZERO) + vr = MACH_EPS * norm * + (ABS (w) + ABS (q) + ABS (x) + ABS (y) + ABS (z)); + { + /* comdiv (x * r - z * ra + q * sa, x * s - z * sa -q * ra, + vr, vi, &h[i][na], &h[i][en]); */ + REAL complex c; + c = (x * r - z * ra + q * sa + I * (x * s - z * sa -q * ra)) / (vr + I * vi); + h[i][na] = creal(c); h[i][en] = cimag(c); + } + if (ABS (x) > ABS (z) + ABS (q)) + { + h[i+1][na] = (-ra - w * h[i][na] + q * h[i][en]) / x; + h[i+1][en] = (-sa - w * h[i][en] - q * h[i][na]) / x; + } + else { + /* comdiv (-r - y * h[i][na], -s - y * h[i][en], z, q, + &h[i+1][na], &h[i+1][en]); */ + REAL complex c; + c = (-r - y * h[i][na] + I * (-s - y * h[i][en])) / (z + I * q); + h[i+1][na] = creal(c); h[i+1][en] = cimag(c); + } + + } /* end wi[i] > 0 */ + } /* end wi[i] >= 0 */ + } /* end i */ + } /* if q < 0 */ + } /* end en */ + + for (i = 0; i < n; i++) /* Eigenvectors for the evalues for */ + if (i < low || i > high) /* rows < low and rows > high */ + for (k = i + 1; k < n; k++) eivec[i][k] = h[i][k]; + + for (j = n - 1; j >= low; j--) + { + m = (j <= high) ? j : high; + if (wi[j] < ZERO) + { + for (l = j - 1, i = low; i <= high; i++) + { + for (y = z = ZERO, k = low; k <= m; k++) + { + y += eivec[i][k] * h[k][l]; + z += eivec[i][k] * h[k][j]; + } + + eivec[i][l] = y; + eivec[i][j] = z; + } + } + else + if (wi[j] == ZERO) + { + for (i = low; i <= high; i++) + { + for (z = ZERO, k = low; k <= m; k++) + z += eivec[i][k] * h[k][j]; + eivec[i][j] = z; + } + } + + } /* end j */ + + return (0); +} + + +static int hqr2 /* compute eigenvalues .......................*/ +/*.IX{hqr2}*/ + (int vec, /* switch for computing evectors*/ + int n, /* Dimension of matrix .........*/ + int low, /* first nonzero row ...........*/ + int high, /* last nonzero row ............*/ + REAL * h[], /* Hessenberg matrix ...........*/ + REAL wr[], /* Real parts of eigenvalues ...*/ + REAL wi[], /* Imaginary parts of evalues ..*/ + REAL * eivec[], /* Matrix of eigenvectors ......*/ + int cnt[] /* Iteration counter ...........*/ + ) +/*====================================================================* + * * + * hqr2 computes the eigenvalues and (if vec != 0) the eigenvectors * + * of an n * n upper Hessenberg matrix. * + * * + *====================================================================* + * * + * Control parameter: * + * ================== * + * vec int vec; * + * = 0 compute eigenvalues only * + * = 1 compute all eigenvalues and eigenvectors * + * * + * Input parameters: * + * ================ * + * n int n; ( n > 0 ) * + * Dimension of mat and eivec, * + * length of the real parts vector wr and of the * + * imaginary parts vector wi of the eigenvalues. * + * low int low; * + * high int high; see balance * + * h REAL *h[n]; * + * upper Hessenberg matrix * + * * + * Output parameters: * + * ================== * + * eivec REAL *eivec[n]; ( bei vec = 1 ) * + * Matrix, which for vec = 1 contains the eigenvectors * + * as follows : * + * For real eigebvalues the corresponding column * + * contains the corresponding eigenvactor, while for * + * complex eigenvalues the corresponding column contains* + * the real part of the eigenvactor with its imaginary * + * part is stored in the subsequent column of eivec. * + * The eigenvactor for the complex conjugate eigenvactor* + * is given by the complex conjugate eigenvactor. * + * wr REAL wr[n]; * + * Real part of the n eigenvalues. * + * wi REAL wi[n]; * + * Imaginary parts of the eigenvalues * + * cnt int cnt[n]; * + * vector of iterations used for each eigenvalue. * + * For a complex conjugate eigenvalue pair the second * + * entry is negative. * + * * + * Return value : * + * ============= * + * = 0 all ok * + * = 4xx Iteration maximum exceeded when computing evalue xx * + * = 99 zero matrix * + * * + *====================================================================* + * * + * functions in use : * + * =================== * + * * + * int hqrvec(): reverse transform for eigenvectors * + * * + *====================================================================* + * * + * Constants used : MACH_EPS, MAXIT * + * ================= * + * * + * Macros : SWAP, ABS, SQRT * + * ========= * + * * + *====================================================================*/ +{ + int i, j; + int na, en, iter, k, l, m; + REAL p = ZERO, q = ZERO, r = ZERO, s, t, w, x, y, z; + + for (i = 0; i < n; i++) + if (i < low || i > high) + { + wr[i] = h[i][i]; + wi[i] = ZERO; + cnt[i] = 0; + } + + en = high; + t = ZERO; + + while (en >= low) + { + iter = 0; + na = en - 1; + + for ( ; ; ) + { + for (l = en; l > low; l--) /* search for small */ + if ( ABS(h[l][l-1]) <= /* subdiagonal element */ + MACH_EPS * (ABS(h[l-1][l-1]) + ABS(h[l][l])) ) break; + + x = h[en][en]; + if (l == en) /* found one evalue */ + { + wr[en] = h[en][en] = x + t; + wi[en] = ZERO; + cnt[en] = iter; + en--; + break; + } + + y = h[na][na]; + w = h[en][na] * h[na][en]; + + if (l == na) /* found two evalues */ + { + p = (y - x) * 0.5; + q = p * p + w; + z = SQRT (ABS (q)); + x = h[en][en] = x + t; + h[na][na] = y + t; + cnt[en] = -iter; + cnt[na] = iter; + if (q >= ZERO) + { /* real eigenvalues */ + z = (p < ZERO) ? (p - z) : (p + z); + wr[na] = x + z; + wr[en] = s = x - w / z; + wi[na] = wi[en] = ZERO; + x = h[en][na]; + r = SQRT (x * x + z * z); + + if (vec) + { + p = x / r; + q = z / r; + for (j = na; j < n; j++) + { + z = h[na][j]; + h[na][j] = q * z + p * h[en][j]; + h[en][j] = q * h[en][j] - p * z; + } + + for (i = 0; i <= en; i++) + { + z = h[i][na]; + h[i][na] = q * z + p * h[i][en]; + h[i][en] = q * h[i][en] - p * z; + } + + for (i = low; i <= high; i++) + { + z = eivec[i][na]; + eivec[i][na] = q * z + p * eivec[i][en]; + eivec[i][en] = q * eivec[i][en] - p * z; + } + } /* end if (vec) */ + } /* end if (q >= ZERO) */ + else /* pair of complex */ + { /* conjugate evalues */ + wr[na] = wr[en] = x + p; + wi[na] = z; + wi[en] = - z; + } + + en -= 2; + break; + } /* end if (l == na) */ + + if (iter >= MAXIT) + { + cnt[en] = MAXIT + 1; + return (en); /* MAXIT Iterations */ + } + + if ( (iter != 0) && (iter % 10 == 0) ) + { + t += x; + for (i = low; i <= en; i++) h[i][i] -= x; + s = ABS (h[en][na]) + ABS (h[na][en-2]); + x = y = (REAL)0.75 * s; + w = - (REAL)0.4375 * s * s; + } + + iter ++; + + for (m = en - 2; m >= l; m--) + { + z = h[m][m]; + r = x - z; + s = y - z; + p = ( r * s - w ) / h[m+1][m] + h[m][m+1]; + q = h[m + 1][m + 1] - z - r - s; + r = h[m + 2][m + 1]; + s = ABS (p) + ABS (q) + ABS (r); + p /= s; + q /= s; + r /= s; + if (m == l) break; + if ( ABS (h[m][m-1]) * (ABS (q) + ABS (r)) <= + MACH_EPS * ABS (p) + * ( ABS (h[m-1][m-1]) + ABS (z) + ABS (h[m+1][m+1])) ) + break; + } + + for (i = m + 2; i <= en; i++) h[i][i-2] = ZERO; + for (i = m + 3; i <= en; i++) h[i][i-3] = ZERO; + + for (k = m; k <= na; k++) + { + if (k != m) /* double QR step, for rows l to en */ + { /* and columns m to en */ + p = h[k][k-1]; + q = h[k+1][k-1]; + r = (k != na) ? h[k+2][k-1] : ZERO; + x = ABS (p) + ABS (q) + ABS (r); + if (x == ZERO) continue; /* next k */ + p /= x; + q /= x; + r /= x; + } + s = SQRT (p * p + q * q + r * r); + if (p < ZERO) s = -s; + + if (k != m) h[k][k-1] = -s * x; + else if (l != m) + h[k][k-1] = -h[k][k-1]; + p += s; + x = p / s; + y = q / s; + z = r / s; + q /= p; + r /= p; + + for (j = k; j < n; j++) /* modify rows */ + { + p = h[k][j] + q * h[k+1][j]; + if (k != na) + { + p += r * h[k+2][j]; + h[k+2][j] -= p * z; + } + h[k+1][j] -= p * y; + h[k][j] -= p * x; + } + + j = (k + 3 < en) ? (k + 3) : en; + for (i = 0; i <= j; i++) /* modify columns */ + { + p = x * h[i][k] + y * h[i][k+1]; + if (k != na) + { + p += z * h[i][k+2]; + h[i][k+2] -= p * r; + } + h[i][k+1] -= p * q; + h[i][k] -= p; + } + + if (vec) /* if eigenvectors are needed ..................*/ + { + for (i = low; i <= high; i++) + { + p = x * eivec[i][k] + y * eivec[i][k+1]; + if (k != na) + { + p += z * eivec[i][k+2]; + eivec[i][k+2] -= p * r; + } + eivec[i][k+1] -= p * q; + eivec[i][k] -= p; + } + } + } /* end k */ + + } /* end for ( ; ;) */ + + } /* while (en >= low) All evalues found */ + + if (vec) /* transform evectors back */ + if (hqrvec (n, low, high, h, wr, wi, eivec)) return (99); + return (0); +} + + +static int norm_1 /* normalize eigenvectors to have one norm 1 .*/ +/*.IX{norm\unt 1}*/ + (int n, /* Dimension of matrix ...........*/ + REAL * v[], /* Matrix with eigenvektors ......*/ + REAL wi[] /* Imaginary parts of evalues ....*/ + ) +/*====================================================================* + * * + * norm_1 normalizes the one norm of the column vectors in v. * + * (special attention to complex vectors in v is given) * + * * + *====================================================================* + * * + * Input parameters: * + * ================ * + * n int n; ( n > 0 ) * + * Dimension of matrix v * + * v REAL *v[]; * + * Matrix of eigenvectors * + * wi REAL wi[]; * + * Imaginary parts of the eigenvalues * + * * + * Output parameter: * + * ================ * + * v REAL *v[]; * + * Matrix with normalized eigenvectors * + * * + * Return value : * + * ============= * + * = 0 all ok * + * = 1 n < 1 * + * * + *====================================================================* + * * + * functions used : * + * ================= * + * REAL comabs(): complex absolute value * + * int comdiv(): complex division * + * * + * Macros : ABS * + * ======== * + * * + *====================================================================*/ +{ + int i, j; + REAL maxi, tr, ti; + + if (n < 1) return (1); + + for (j = 0; j < n; j++) + { + if (wi[j] == ZERO) + { + maxi = v[0][j]; + for (i = 1; i < n; i++) + if (ABS (v[i][j]) > ABS (maxi)) maxi = v[i][j]; + + if (maxi != ZERO) + { + maxi = ONE / maxi; + for (i = 0; i < n; i++) v[i][j] *= maxi; + } + } + else + { + tr = v[0][j]; + ti = v[0][j+1]; + for (i = 1; i < n; i++) + /* if ( comabs (v[i][j], v[i][j+1]) > comabs (tr, ti) ) */ + if (cabs(v[i][j] + I * v[i][j+1]) > cabs(tr + I * ti)) + { + tr = v[i][j]; + ti = v[i][j+1]; + } + + if (tr != ZERO || ti != ZERO) + for (i = 0; i < n; i++) { + /* comdiv (v[i][j], v[i][j+1], tr, ti, &v[i][j], &v[i][j+1]); */ + REAL complex c; + c = (v[i][j] + I * v[i][j+1]) / (tr + I * ti); + v[i][j] = creal(c); v[i][j+1] = cimag(c); + } + + j++; /* raise j by two */ + } + } + return (0); +} + +/*.BA*/ + +static int eigen /* Compute all evalues/evectors of a matrix ..*/ +/*.IX{eigen}*/ + ( + int vec, /* switch for computing evectors ...*/ + int ortho, /* orthogonal Hessenberg reduction? */ + int ev_norm, /* normalize Eigenvectors? .........*/ + int n, /* size of matrix ..................*/ + REAL * mat[], /* input matrix ....................*/ + REAL * eivec[], /* Eigenvectors ....................*/ + REAL valre[], /* real parts of eigenvalues .......*/ + REAL valim[], /* imaginary parts of eigenvalues ..*/ + int cnt[] /* Iteration counter ...............*/ + ) +/*====================================================================* + * * + * The function eigen determines all eigenvalues and (if desired) * + * all eigenvectors of a real square n * n matrix via the QR method* + * in the version of Martin, Parlett, Peters, Reinsch and Wilkinson.* + * * + *====================================================================* +.BE*) + * * + * Literature: * + * =========== * + * 1) Peters, Wilkinson: Eigenvectors of real and complex * + * matrices by LR and QR triangularisations, * + * Num. Math. 16, p.184-204, (1970); [PETE70]; contribution * + * II/15, p. 372 - 395 in [WILK71]. * + * 2) Martin, Wilkinson: Similarity reductions of a general * + * matrix to Hessenberg form, Num. Math. 12, p. 349-368,(1968)* + * [MART 68]; contribution II,13, p. 339 - 358 in [WILK71]. * + * 3) Parlett, Reinsch: Balancing a matrix for calculations of * + * eigenvalues and eigenvectors, Num. Math. 13, p. 293-304, * + * (1969); [PARL69]; contribution II/11, p.315 - 326 in * + * [WILK71]. * + * * + *====================================================================* + * * + * Control parameters: * + * =================== * + * vec int vec; * + * call for eigen : * + * = 0 compute eigenvalues only * + * = 1 compute all eigenvalues and eigenvectors * + * ortho flag that shows if transformation of mat to * + * Hessenberg form shall be done orthogonally by * + * `orthes' (flag set) or elementarily by `elmhes' * + * (flag cleared). The Householder matrices used in * + * orthogonal transformation have the advantage of * + * preserving the symmetry of input matrices. * + * ev_norm flag that shows if Eigenvectors shall be * + * normalized (flag set) or not (flag cleared) * + * * + * Input parameters: * + * ================ * + * n int n; ( n > 0 ) * + * size of matrix, number of eigenvalues * + * mat REAL *mat[n]; * + * matrix * + * * + * Output parameters: * + * ================== * + * eivec REAL *eivec[n]; ( bei vec = 1 ) * + * matrix, if vec = 1 this holds the eigenvectors * + * thus : * + * If the jth eigenvalue of the matrix is real then the * + * jth column is the corresponding real eigenvector; * + * if the jth eigenvalue is complex then the jth column * + * of eivec contains the real part of the eigenvector * + * while its imaginary part is in column j+1. * + * (the j+1st eigenvector is the complex conjugate * + * vector.) * + * valre REAL valre[n]; * + * Real parts of the eigenvalues. * + * valim REAL valim[n]; * + * Imaginary parts of the eigenvalues * + * cnt int cnt[n]; * + * vector containing the number of iterations for each * + * eigenvalue. (for a complex conjugate pair the second * + * entry is negative.) * + * * + * Return value : * + * ============= * + * = 0 all ok * + * = 1 n < 1 or other invalid input parameter * + * = 2 insufficient memory * + * = 10x error x from balance() * + * = 20x error x from elmh() * + * = 30x error x from elmtrans() (for vec = 1 only) * + * = 4xx error xx from hqr2() * + * = 50x error x from balback() (for vec = 1 only) * + * = 60x error x from norm_1() (for vec = 1 only) * + * * + *====================================================================* + * * + * Functions in use : * + * =================== * + * * + * static int balance (): Balancing of an n x n matrix * + * static int elmh (): Transformation to upper Hessenberg form * + * static int elmtrans(): intialize eigenvectors * + * static int hqr2 (): compute eigenvalues/eigenvectors * + * static int balback (): Reverse balancing to obtain eigenvectors * + * static int norm_1 (): Normalize eigenvectors * + * * + * void *vmalloc(): allocate vector or matrix * + * void vmfree(): free list of vectors and matrices * + * * + *====================================================================* + * * + * Constants used : NULL, BASIS * + * =================== * + * * + *====================================================================*/ +{ + int i; + int low, high, rc; + REAL *scale, + *d = NULL; + void *vmblock; + + if (n < 1) return (1); /* n >= 1 ............*/ + + if (valre == NULL || valim == NULL || mat == NULL || cnt == NULL) + return (1); + + for (i = 0; i < n; i++) + if (mat[i] == NULL) return (1); + + for (i = 0; i < n; i++) cnt[i] = 0; + + if (n == 1) /* n = 1 .............*/ + { + eivec[0][0] = ONE; + valre[0] = mat[0][0]; + valim[0] = ZERO; + return (0); + } + + if (vec) + { + if (eivec == NULL) return (1); + for (i = 0; i < n; i++) + if (eivec[i] == NULL) return (1); + } + + vmblock = vminit(); + scale = (REAL *)vmalloc(vmblock, VEKTOR, n, 0); + if (! vmcomplete(vmblock)) /* memory error */ + return 2; + + if (vec && ortho) /* with Eigenvectors */ + { /* and orthogonal */ + /* Hessenberg reduction? */ + d = (REAL *)vmalloc(vmblock, VEKTOR, n, 0); + if (! vmcomplete(vmblock)) + { + vmfree(vmblock); + return 1; + } + } + + /* balance mat for nearly */ + rc = balance (n, mat, scale, /* equal row and column */ + &low, &high, BASIS); /* one norms */ + if (rc) + { + vmfree(vmblock); + return (100 + rc); + } + + if (ortho) + rc = orthes(n, low, high, mat, d); + else + rc = elmhes (n, low, high, mat, cnt); /* reduce mat to upper */ + if (rc) /* Hessenberg form */ + { + vmfree(vmblock); + return (200 + rc); + } + + if (vec) /* initialize eivec */ + { + if (ortho) + rc = orttrans(n, low, high, mat, d, eivec); + else + rc = elmtrans (n, low, high, mat, cnt, eivec); + if (rc) + { + vmfree(vmblock); + return (300 + rc); + } + } + + rc = hqr2 (vec, n, low, high, mat, /* execute Francis QR */ + valre, valim, eivec, cnt); /* algorithm to obtain */ + if (rc) /* eigenvalues */ + { + vmfree(vmblock); + return (400 + rc); + } + + if (vec) + { + rc = balback (n, low, high, /* reverse balancing if */ + scale, eivec); /* eigenvaectors are to */ + if (rc) /* be determined */ + { + vmfree(vmblock); + return (500 + rc); + } + if (ev_norm) + rc = norm_1 (n, eivec, valim); /* normalize eigenvectors */ + if (rc) + { + vmfree(vmblock); + return (600 + rc); + } + } + + vmfree(vmblock); /* free buffers */ + + + return (0); +} + +/* -------------------------- END feigen.c -------------------------- */ + +/* _a[0..n^2-1] is a real general matrix. On return, evalr store the + real part of eigenvalues and evali the imgainary part. If _evec is + not a NULL pointer, eigenvectors will be stored there. */ +int n_eigeng(double *_a, int n, double *evalr, double *evali, double *_evec) +{ + double **a, **evec = 0; + int i, j, *cnt; + a = (double**)calloc(n, sizeof(void*)); + if (_evec) evec = (double**)calloc(n, sizeof(void*)); + cnt = (int*)calloc(n, sizeof(int)); + for (i = 0; i < n; ++i) { + a[i] = _a + i * n; + if (_evec) evec[i] = _evec + i * n; + } + eigen(_evec? 1 : 0, 0, 1, n, a, evec, evalr, evali, cnt); + if (_evec) { + double tmp; + for (j = 0; j < n; ++j) { + tmp = 0.0; + for (i = 0; i < n; ++i) tmp += SQR(evec[i][j]); + tmp = SQRT(tmp); + for (i = 0; i < n; ++i) evec[i][j] /= tmp; + } + } + free(a); free(evec); free(cnt); + return 0; +} + +#ifdef LH3_EIGENG_MAIN +int main(void) +{ + static double a[5][5] = {{1.0, 6.0, -3.0, -1.0, 7.0}, + {8.0, -15.0, 18.0, 5.0, 4.0}, + {-2.0, 11.0, 9.0, 15.0, 20.0}, + {-13.0, 2.0, 21.0, 30.0, -6.0}, + {17.0, 22.0, -5.0, 3.0, 6.0}}; + + static double b[5][5]={ {10.0,1.0,2.0,3.0,4.0}, + {1.0,9.0,-1.0,2.0,-3.0}, + {2.0,-1.0,7.0,3.0,-5.0}, + {3.0,2.0,3.0,12.0,-1.0}, + {4.0,-3.0,-5.0,-1.0,15.0}}; + + double *mem, *evalr, *evali, *evec; + int i, j; + mem = (double*)calloc(5 * 5 + 10, sizeof(double)); + evec = mem; + evalr = evec + 25; + evali = evalr + 5; +// n_eigeng(a[0], 5, evalr, evali, evec); + n_eigeng(a[0], 5, evalr, evali, 0); + for (i = 0; i < 5; ++i) { + printf("%le + %le J\n", evalr[i], evali[i]); + } + + printf("\n\n"); + n_eigeng(b[0], 5, evalr, evali, evec); + for (i = 0; i <= 4; i++) + printf("%13.7e + %13.7e J\n", evalr[i], evali[i]); + printf("\n"); + for (i = 0; i <= 4; i++) { + for (j = 0; j <= 4; j++) + printf("%12.6e ", evec[i*5 + j]); + printf("\n"); + } + printf("\n"); + + free(mem); + return 0; +} +#endif diff --git a/libs/DISTlib/source/file_reader.c b/libs/DISTlib/source/file_reader.c new file mode 100644 index 000000000..21966e6e7 --- /dev/null +++ b/libs/DISTlib/source/file_reader.c @@ -0,0 +1,526 @@ +#include +#include +#include +#include +#include +#include "helper.h" +#include "distinterface.h" +#include "outputdist.h" +#include "file_reader.h" +#include "distgeneration.h" +//#include "round_near.c" + + +int readfile(const char* filename){ + dist->incoordtype=2; + FILE *file = fopen ( filename, "r" ); + double** table = malloc(MAX_ROWS * sizeof(double*)); // allocate the rows + for (int r = 0; r < MAX_ROWS; ++r){ + table[r] = malloc(MAX_COLUMNS * sizeof(double)); // allocate the columns + } + + int linecount = -1; + int numcolum, index; + char columns[MAX_COLUMNS][MAX_LENGTH]; + char units[MAX_COLUMNS][MAX_LENGTH]; + char line [ MAX_LENGTH*50 ]; /* or other suitable maximum line size */ + char tosplit [ MAX_LENGTH*10 ]; + char value_s[MAX_LENGTH], shorty[MAX_LENGTH]; + char * parameter = (char*)malloc(MAX_LENGTH*sizeof(char)); + char * parameter_tmp = (char*)malloc(MAX_LENGTH*sizeof(char)); + char * unit =(char*)malloc(MAX_LENGTH*sizeof(char)); + char * unit_tmp =(char*)malloc(MAX_LENGTH*sizeof(char)); + char * at =(char*)malloc(MAX_LENGTH*sizeof(char)); + double value; + char *e ; + if ( file != NULL ){ + while ( fgets ( line, sizeof line, file ) != NULL ){ /* read a line */ + double multifactor = 1; + int k =0; + if(strncmp(line, "@", 1)==0){ + while( line[k] ) { + line[k] = (tolower(line[k])); + k++; + } + sscanf( line, "%s %s %s", at, parameter_tmp, value_s); + if(strncmp(value_s, "%", 1)!=0){ + unit_tmp =strstr(parameter_tmp, "["); + + if(unit_tmp != NULL){ + strcpy(unit,unit_tmp); + e = strchr(parameter_tmp, '['); + index = (int)(e - parameter_tmp); + strncpy(parameter, parameter_tmp, index); + } + else{ + strcpy(unit, "nounit"); + strcpy(parameter, parameter_tmp); + } + + value = strtod(value_s,NULL); + + if(strcmp(parameter, "energy0")==0){ + multifactor = getEnergyUnit(unit); + dist->ref->e0=value*multifactor; + } + else if(strcmp(parameter, "pc0")==0){ + multifactor = getEnergyUnit(unit); + dist->ref->pc0=value*multifactor; + } + else if(strcmp(parameter, "a0")==0){ + dist->ref->a0=round(value); + } + else if(strcmp(parameter, "z0")==0){ + dist->ref->z0=round(value); + } + else if(strcmp(parameter, "mass0")==0){ + multifactor = getEnergyUnit(unit); + dist->ref->mass0=value*multifactor; + } + else if(strcmp(parameter, "charge0")==0){ + dist->ref->charge0=round(value); + } + else{ + issue_error2("Not recogniced parameter!", parameter); + } + } + } + else if(strncmp(line, "!", 1)==0 || strncmp(line, "$", 1)==0){ + printf("Comment line %s \n", line ); + } + else if(strncmp(line, "!", 1)!=0 && linecount==-1){ + + strcpy(tosplit, line); + numcolum = splitline(tosplit, columns, units); + linecount++; + } + else if(linecount>=0 && strncmp(line, "$", 1)!=0){ + add2table(table, line, linecount); + linecount++; + } + else{ + issue_error2("Something went wrong with reading line", line); + } + + } + fclose ( file ); + } + else + { + perror ( filename ); + } + allocateincoord(linecount); //alocates the memory + fillcoordstructure(numcolum,columns, units, table); //fils the coordinate structures + calculaterefparam(); //Sets up the necessary ref parameters and checks for consistensy. + convert2standard(); // converts to the units used internal , x, px, y, py, zeta, deltap + checkinputtype(); + free_readin_memory(table); + return 0; +} +void checkinputtype(){ + if(dist->ref->typeused[0]==dist->ref->typeused[1]==dist->ref->typeused[2]==dist->ref->typeused[3]==dist->ref->typeused[4]==dist->ref->typeused[5]) + dist->incoordtype=dist->ref->typeused[0]; + else if(dist->ref->typeused[4]==2==dist->ref->typeused[5]==2 && dist->ref->typeused[0]==dist->ref->typeused[1]==dist->ref->typeused[2]==dist->ref->typeused[3]==0) + dist->incoordtype=3; //Mixed with action angle coordinates + else if(dist->ref->typeused[4]==2==dist->ref->typeused[5]==2 && dist->ref->typeused[0]==dist->ref->typeused[1]==dist->ref->typeused[2]==dist->ref->typeused[3]==1) + dist->incoordtype=4; //Mixed with normalised coordinates + else + issue_error("Non comptable types have been entered."); + + if(dist->incoordtype==0){ + + } + else if(dist->incoordtype==1){ + + } + else if(dist->incoordtype==3){ + + } + else if(dist->incoordtype==4){ + + } + +} + +// This checks that reference energy is only set once and sets up the mass in case it is not a column +void calculaterefparam(){ + if(dist->ref->mass0 == 0){ + issue_error("A mass is needed! \n"); + } + if(dist->ref->pc0 == 0 && dist->ref->e0==0){ + issue_error("A energy is needed! \n"); + } + if(dist->ref->pc0 > 0 && dist->ref->e0 > 0){ + issue_error("Can't set both energy and momentum! \n"); + } + if(dist->ref->pc0 > 0){ + dist->ref->e0 = momentum2energy(dist->ref->pc0,dist->ref->mass0); + } + if(dist->ref->e0 > 0){ + dist->ref->pc0 = energy2momentum(dist->ref->e0,dist->ref->mass0); + } + + dist->ref->beta0 = (dist->ref->pc0)/(dist->ref->e0); + + if(dist->incoord[0]->mass < 1e-16){ //enough to check the first + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->mass = dist->ref->mass0; //Gives all the same mass. + } + } + if(dist->incoord[0]->a ==0){ //enough to check the first + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->a = dist->ref->a0; //Gives all the same mass. + } + } + if(abs(dist->incoord[0]->z) ==0){ //enough to check the first + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->z = dist->ref->z0; //Gives all the same mass. + } + } +} + +//Converts to the internal used canonical +void convert2standard(){ + for(int i=0; itotincoord; i++){ + for(int j=0; j<6; j++){ + dist->incoord[i]->coord[j]= dist->incoord[i]->readin[j]; // have to have checks here..! + } + } + + if(dist->ref->en_like==-1){ + issue_info("No energy variable is set. Assume 0 deviation from reference energy \n"); + } + else if(dist->ref->en_like==0){ //energy + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->coord[5] = (momentum2energy(dist->incoord[i]->readin[6], dist->incoord[i]->mass)-(dist->ref->pc0))/(dist->ref->pc0); + } + } + else if(dist->ref->en_like==1){ //momentum + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->coord[5] = ((dist->incoord[i]->readin[7], dist->incoord[i]->mass)-(dist->ref->pc0))/(dist->ref->pc0); + } + } + else if(dist->ref->en_like==2){ //psigma + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->coord[5] = psigma2deltap(dist->incoord[i]->readin[8], dist->ref->beta0); + } + } + else if(dist->ref->en_like==3){ //pt + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->coord[5] = pt2deltap(dist->incoord[i]->readin[9], dist->ref->beta0); + } + } + + if(dist->ref->time_like==-1){ + issue_info("No time variable is set. Assume 0 deviation \n"); + } + else if(dist->ref->time_like==2){ //sigma + double beta, pc; + for(int i=0; i< dist->totincoord; i++){ + pc = dist->ref->pc0+dist->incoord[i]->coord[5]; + beta = (pc)/momentum2energy(pc, dist->incoord[i]->mass); + dist->incoord[i]->coord[4] = sigma2zeta(dist->incoord[i]->readin[10], dist->ref->beta0, beta); + } + } + else if(dist->ref->time_like==3){ //t or tau + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->coord[4] = tau2zeta(dist->incoord[i]->readin[11], dist->ref->beta0); + } + } + if(dist->ref->ang_like==-1){ + issue_error("No angle variable is set. \n"); + } + else if(dist->ref->ang_like==1){ //t or tau + for(int i=0; i< dist->totincoord; i++){ + dist->incoord[i]->coord[1] = dist->incoord[i]->coord[6]/(1+dist->incoord[i]->coord[5]); + dist->incoord[i]->coord[3] = dist->incoord[i]->coord[7]/(1+dist->incoord[i]->coord[5]); + } + } +} + +void allocateincoord(int linecount){ + dist->incoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); + dist->outcoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); + dist->totincoord = linecount; + for(int i=0; i< dim; i++){ + dist->ref->typeused[i]=-1; + } + for(int i=0; iincoord[i] = (struct coordinates*)malloc(sizeof(struct coordinates)); + + dist->incoord[i]->coord = (double*)malloc(dim*sizeof(double)); + dist->incoord[i]->readin = (double*)malloc(32*sizeof(double)); + + dist->incoord[i]->mass=0; + dist->incoord[i]->a=0; + dist->incoord[i]->z=0; + + dist->outcoord[i] = (struct coordinates*)malloc(sizeof(struct coordinates)); + dist->outcoord[i]->coord = (double*)malloc(dim*sizeof(double)); + + } + dist->isallocated =1; +} +void deallocateincoord(){ + + + for(int i=0; itotincoord; i++){ + + // free(dist->incoord[i]->coord); + // free(dist->incoord[i]->readin); + // free(dist->incoord[i]); + // free(dist->outcoord[i]->coord); +// free(dist->outcoord[i]); + + } + + free(dist->incoord); + free(dist->outcoord); + //dist->ref->typeused = (int*)malloc(dim*sizeof(int)); + dist->isallocated =0; + dist->totincoord = -1; +} + + + +void setreadin(int coordorder, int column, double ** table, double multifactor, int realorder){ + if(dist->ref->typeused[realorder]!=-1){ + issue_error("Only allowed to set one coordinates index ones"); + } + + if(coordorder < 20){ + dist->ref->typeused[realorder] = 2; + } + else if(coordorder >= 20 <=25){ + dist->ref->typeused[realorder] = 1; + } + else{ + dist->ref->typeused[realorder] = 0; + } + for(int i=0;i < dist->totincoord; i++){ + dist->incoord[i]->readin[coordorder] = multifactor*table[i][column]; + } +} + + +int splitline(char* line, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH] ){ + char* word; + int count = 0; + int k=0; + + while( line[k] ) { + line[k] = (tolower(line[k])); + k++; + } + + word = strtok(line, " "); + if(strncmp("*", word, 1)!=0){ + add2internaltab(word, columns, units, count); + count++; + } + while ((word = strtok(NULL, " ")) != NULL){ + add2internaltab(word, columns, units, count); + count ++; + } + + return count; +} +void add2internaltab(char* word, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH] , int count){ + + + char * unit =(char*)malloc(MAX_LENGTH*sizeof(char)); + char * word_tmp =(char*)malloc(MAX_LENGTH*sizeof(char)); + char * columnname =(char*)malloc(MAX_LENGTH*sizeof(char)); + char *e; + int index; + strcpy(word_tmp, word); + unit =strstr(word_tmp, "["); + if(unit != NULL){ + + strcpy(units[count], unit); + e = strchr(word_tmp, '['); + index = (int)(e - word_tmp); + strncpy(columns[count], word_tmp, index); + columns[count][index]='\0'; + } + else{ + strcpy(units[count], "nounit"); + strcpy(columns[count], word_tmp); + } +} + +void add2table(double ** table, char* line, int linenum){ + char* word; + int count = 0; + word = strtok(line, " "); + table[linenum][0] = atof(word); + count++; + + /* the following loop gets the rest of the words until the + * end of the message */ + while ((word = strtok(NULL, " ")) != NULL){ + table[linenum][count] = atof(word); + count ++; + } +} + +void fillcoordstructure(int numcolum, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH], double **table){ + + double multifactor; + for(int i=0; i< numcolum; i++){ + if(strcmp(columns[i], "x")==0){ + multifactor = getMetricUnit(units[i]); + setreadin(0, i, table, multifactor, 0); + } + else if(strcmp(columns[i], "px")==0){ + checkifangset(0); + setreadin(1, i, table,1, 1); + } + else if(strcmp(columns[i], "y")==0){ + multifactor = getMetricUnit(units[i]); + setreadin(2, i, table, multifactor, 2); + } + else if(strcmp(columns[i], "py")==0){ + checkifangset(0); + setreadin(3, i, table,1, 3); + } + else if(strcmp(columns[i], "zeta")==0){ + checkiftimeset(5); + multifactor = getMetricUnit(units[i]); + setreadin(4, i, table, multifactor, 4); + } + else if(strcmpnl(columns[i], "deltap")==0){ + checkifenergyset(5); + setreadin(5, i, table, 1, 5); + } + else if(strcmp(columns[i], "energy")==0){ + multifactor = getEnergyUnit(units[i]); + checkifenergyset(0); + //setnonstandard(0, i, table); + setreadin(6, i, table, 1, 5); + } + else if(strcmp(columns[i], "pc")==0){ + multifactor = getEnergyUnit(units[i]); + checkifenergyset(1); + //setnonstandard(1, i, table); + setreadin(7, i, table, 1, 5); + } + else if(strcmp(columns[i], "psigma")==0){ + checkifenergyset(2); + // setnonstandard(2, i, table); + setreadin(8, i, table, 1, 5); + } + else if(strcmp(columns[i], "ptau")==0 || strcmp(columns[i], "pt")==0 ){ + checkifenergyset(3); + // setnonstandard(3, i, table); + setreadin(9, i, table, 1, 4); + } + else if(strcmp(columns[i], "sigma")==0){ + checkiftimeset(2); + multifactor = getMetricUnit(units[i]); + //setnonstandard(4, i, table); + setreadin(10, i, table, 1, 4); + } + else if(strcmp(columns[i], "tau")==0 || strcmp(columns[i], "t")==0 ){ + checkiftimeset(3); + multifactor = getMetricUnit(units[i]); + //setnonstandard(5, i, table); + setreadin(11, i, table, 1, 4); + } + else if(strcmp(columns[i], "xp")==0){ + checkifangset(1); + //setnonstandard(6, i, table); + setreadin(12, i, table, 1, 1); + } + else if(strcmp(columns[i], "yp")==0){ + checkifangset(1); + //setnonstandard(7, i, table); + setreadin(13, i, table, 1, 3); + } + else if(strcmp(columns[i], "jx")==0) + setreadin(20, i, table, 1, 0); + else if(strcmp(columns[i], "phix")==0) + setreadin(21, i, table, 1, 1); + else if(strcmp(columns[i], "jy")==0) + setreadin(22, i, table, 1, 2); + else if(strcmp(columns[i], "phiy")==0) + setreadin(23, i, table, 1, 3); + else if(strcmp(columns[i], "jz")==0) + setreadin(24, i, table, 1, 4); + else if(strcmp(columns[i], "phiz")==0) + setreadin(25, i, table, 1, 5); + else if(strcmp(columns[i], "xn")==0) + setreadin(26, i, table, 1, 0); + else if(strcmp(columns[i], "pxn")==0) + setreadin(27, i, table, 1, 1); + else if(strcmp(columns[i], "yn")==0) + setreadin(28, i, table, 1, 2); + else if(strcmp(columns[i], "pyn")==0) + setreadin(29, i, table, 1, 3); + else if(strcmp(columns[i], "zn")==0) + setreadin(30, i, table, 1, 4); + else if(strcmp(columns[i], "pzn")==0) + setreadin(31, i, table, 1, 5); + else if(strcmp(columns[i], "mass")==0){ + multifactor = getEnergyUnit(units[i]); + for(int j=0;j < dist->totincoord; j++){ + dist->incoord[j]->mass = multifactor*table[j][i]; + } + } + else if(strcmp(columns[i], "a")==0){ + for(int j=0;j < dist->totincoord; j++){ + dist->incoord[j]->a = table[j][i]; + } + } + else if(strcmp(columns[i], "z")==0){ + for(int j=0;j < dist->totincoord; j++){ + dist->incoord[j]->z = table[j][i]; + } + } + } +} +double getEnergyUnit(char *unit){ + double multif; + if(strcmp(unit, "[ev]")==0){ + multif = 1e-9; + } + else if(strcmp(unit, "[kev]")==0){ + multif = 1e-6; + } + else if(strcmp(unit, "[mev]")==0){ + multif = 1e-3; + } + else if(strcmp(unit, "[gev]")==0 || strcmp(unit, "nounit")==0){ + multif = 1; + } + else if(strcmp(unit, "[tev]")==0){ + multif = 1e3; + } + else{ + issue_error("Unknow unit. Must be of the form [Mev]"); + } + return multif; +} + +double getMetricUnit(char *unit){ + double multif; + + if(strcmp(unit, "[mm]")==0){ + multif = 1e-3; + } + else if(strcmp(unit, "[m]")==0 || strcmp(unit, "nounit")==0 ){ + multif = 1; + } + else{ + issue_warning("Unknow unit. Must be of the form [Mev]"); + } + return multif; +} + +void free_readin_memory(double **table){ + + for (int r = 0; r < MAX_ROWS; ++r){ + free(table[r]); + } + free(table); +} \ No newline at end of file diff --git a/libs/DISTlib/source/file_reader.h b/libs/DISTlib/source/file_reader.h new file mode 100644 index 000000000..de5cda18c --- /dev/null +++ b/libs/DISTlib/source/file_reader.h @@ -0,0 +1,18 @@ +int readfile(const char* filenamet); +int splitline(char* line, char columns[MAX_COLUMNS][MAX_LENGTH],char units[MAX_COLUMNS][MAX_LENGTH] ); +void add2table(double ** table, char* line, int linenum ); +void allocateincoord(int linecount); +void setreadin(int coordorder, int column, double ** table, double multifactor, int realcoord); +void checkifenergyset(int entype); +void checkiftimeset(int entype); +void checkifangset(int entype); +void calculaterefparam(); +void convert2standard(); +void add2internaltab(char* word, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH], int count) ; +double getEnergyUnit(char *unit); +double getMetricUnit(char *unit); +void fillcoordstructure(int numcolum, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH], double **table); +void free_readin_memory(double **table); +void checkinputtype(); +void deallocateincoord(); +void allocateincoord(); \ No newline at end of file diff --git a/libs/DISTlib/source/helper.c b/libs/DISTlib/source/helper.c new file mode 100644 index 000000000..6b28b653e --- /dev/null +++ b/libs/DISTlib/source/helper.c @@ -0,0 +1,148 @@ +#include +#include +#include +#include +#include +#include "helper.h" +#include "distinterface.h" +#include "distgeneration.h" +#include "outputdist.h" + +/*This function converts from canoncial to sixtrack tracking variables*/ +void canonical2six(double *canonical, double beta0, double pc0, double mass0, double mass, double *coord){ + double deltap = canonical[5]; + double beta = (pc0+deltap)/momentum2energy((pc0+deltap), mass); + double rv = beta0/beta; + double factor = 1e3; + coord[0] = canonical[0]*factor; + coord[1] = canonical[1]*(factor/(1+deltap)); + coord[2] = canonical[2]*factor; + coord[3] = canonical[3]*(factor/(1+deltap)); + coord[4] = canonical[4]*(factor*rv); + coord[5] = canonical[5]; +} +/* Compare two strings s1 and s2, assuming either is + * terminated by \n or a NULL, A match + * returns 0, a non-match returns 1. + */ +int strcmpnl(const char *s1, const char *s2) +{ + char s1c; + char s2c; + do + { + s1c = *(s1++); + s2c = *(s2++); + if (s1c == '\n') + s1c = 0; + if (s2c == '\n') + s2c = 0; + if (s1c != s2c) + return 1; + } + while (s1c); /* already checked *s2 is equal */ + return 0; +} + +void issue_error(const char* t1){ + printf("+=+=+= fatal: %s\n",t1); + exit(1); +} + +void issue_error2(const char* t1,const char* t2){ + printf("+=+=+= fatal: %s %s\n",t1, t2); + exit(1); +} + +void issue_warning(const char* t1){ + printf( "+=+=+= warning: %s\n",t1); +} + +void issue_info(const char* t1){ + printf( "Info: %s\n",t1); +} + +// This function prevents several time variables from beeing set +void checkiftimeset(int entype){ + if(dist->ref->time_like==-1) + dist->ref->time_like=entype; + else { + issue_error("Only allowed 1 type of time variable!"); + } +} + +void checkifangset(int entype){ + if(dist->ref->ang_like==-1 || dist->ref->ang_like==entype ) + dist->ref->ang_like=entype; + else { + issue_error("Only allowed 1 type of angle variable!"); + } +} +// This function prevents several energy from beeing set +void checkifenergyset(int entype){ + if(dist->ref->en_like==-1) + dist->ref->en_like=entype; + else { + issue_error("Only allowed 1 type of energy variable!"); + } +} + +double psigma2deltap(double psigma, double beta0 ){ + return (sqrt(pow(psigma*beta0,2) +2*psigma +1)-1); +} + +double sigma2zeta(double sigma, double beta0, double beta){ + return sigma*(beta/beta0); +} + +double tau2zeta(double tau, double beta){ + return tau*beta; +} + +double pt2deltap(double pt, double beta0 ){ + return (sqrt(pow(pt,2) +2*pt*beta0 +1)-1); +} + +double momentum2energy(double momentum, double mass){ + return sqrt(pow(momentum,2)+pow(mass,2)); +} + +double energy2momentum(double energy, double mass){ + double momsq = pow(energy,2)-pow(mass,2); + if(momsq <= 0){ + issue_error("Energy must be larger than mass!"); + } + return sqrt(momsq); +} + +void solve2by2eq(double a1, double b1, double c1, double a2, double b2, double c2, double *x){ + double det = a1*b2-b1*a2; + double dx = c1*b2-b1*c2; + double dy = a1*c2-c1*a2; + x[0] = (dx/det); + x[1] = (dy/det); + +} +void mtrx_vector_mult_pointer(int mp, int np, double **mtrx_a, double mtrx_b[6], double result[6]) +{ + + int m = mp; + int n = np; + result[0]=0; + result[1]=0; + result[2]=0; + result[3]=0; + result[4]=0; + result[5]=0; + + register int i=0, j=0, k=0; + for (i = 0; i < mp; i++) + { + for (k = 0; k < np; k++) + { + result [i] += mtrx_a [i][k] * mtrx_b [k]; + + } + } + +} \ No newline at end of file diff --git a/libs/DISTlib/source/helper.h b/libs/DISTlib/source/helper.h new file mode 100644 index 000000000..bb189d4b2 --- /dev/null +++ b/libs/DISTlib/source/helper.h @@ -0,0 +1,17 @@ +#define MAX_COLUMNS 100 +#define MAX_LENGTH 100 +#define MAX_ROWS 100000 +void canonical2six(double *canonical, double beta0, double pc0, double mass0, double mass, double *coord); +double momentum2energy(double momentum, double mass); +double energy2momentum(double energy, double mass); +double pt2deltap(double pt, double beta0 ); +double psigma2deltap(double psigma, double beta0 ); +void issue_error(const char* t1); +void issue_error2(const char* t1,const char* t2); +void issue_warning(const char* t1); +void issue_info(const char* t1); +int strcmpnl (const char *s1, const char *s2); +double sigma2zeta(double sigma, double beta0, double beta); +double tau2zeta(double tau, double beta); +void solve2by2eq(double a1, double b1, double c1, double a2, double b2, double c2, double *x); +void mtrx_vector_mult_pointer(int mp, int np, double **mtrx_a, double mtrx_b[6], double result[6]); \ No newline at end of file diff --git a/libs/DISTlib/source/outputdist.c b/libs/DISTlib/source/outputdist.c new file mode 100644 index 000000000..27be1adbd --- /dev/null +++ b/libs/DISTlib/source/outputdist.c @@ -0,0 +1,33 @@ +#include +#include +#include +#include +#include "helper.h" +#include "distinterface.h" +#include "distgeneration.h" +#include "outputdist.h" + + +void print2file(const char* nameoffile){ + + FILE * fp; + /* open the file for writing*/ + fp = fopen (nameoffile,"w"); + fprintf (fp, "@ mass0 %f \n",dist->ref->mass0); + fprintf (fp, "@ charge0 %d \n",dist->ref->charge0); + fprintf (fp, "@ z0 %d \n",dist->ref->z0); + fprintf (fp, "@ a0 %d \n",dist->ref->a0); + fprintf (fp, "@ pc0 %f \n",dist->ref->pc0); + + fprintf(fp, "x px y py zeta deltap \n"); + if(dist->isDistrcalculated ==0){ + gensixcanonical(); + } + for(int i=0; itotoutcoord; i++){ + fprintf(fp, "%.9e %.9e %.9e %.9e %.9e %.9e \n", dist->outcoord[i]->coord[0],dist->outcoord[i]->coord[1],dist->outcoord[i]->coord[2], + dist->outcoord[i]->coord[3],dist->outcoord[i]->coord[4],dist->outcoord[i]->coord[5]); + } + + /* close the file*/ + fclose (fp); +} \ No newline at end of file diff --git a/libs/DISTlib/source/outputdist.h b/libs/DISTlib/source/outputdist.h new file mode 100644 index 000000000..393c89061 --- /dev/null +++ b/libs/DISTlib/source/outputdist.h @@ -0,0 +1,10 @@ +//void print2file(); +//void printdistsettings_(int *ndist); +//int getnumberdist_(); +//int getupperbound(); +//void getcoordvectors_(double *x, double *xp, double *y, double *yp, double *sigma, double *delta); +//void getcoord_(double coordinate[6], int initial ); +//void getcoord_normalized(double coordinate_normalized[6], int initial ); +//void getcoord2_(double coordinate[6], int *initial ); +//void print2fort13(); +void print2file(const char* nameoffile); \ No newline at end of file From 85bebea0a8b4f5f368314799112a3511c9ab9f59 Mon Sep 17 00:00:00 2001 From: Tobias Date: Mon, 23 Mar 2020 18:32:48 +0100 Subject: [PATCH 03/20] This seems to fix the compilation. --- Makefile_c | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/Makefile_c b/Makefile_c index b38f2e772..b0d353ae0 100644 --- a/Makefile_c +++ b/Makefile_c @@ -16,16 +16,19 @@ # | $Id$ # | -vpath %.h src -vpath %.c src +vpath %.h src libs/DISTlib/source +vpath %.c src libs/DISTlib/source vpath %.d $(OBJDIR) -CC_HDR := $(notdir $(wildcard src/mad_*.h)) -CC_SRC := $(notdir $(wildcard src/mad_*.c)) +CC_HDR :=$(notdir $(wildcard libs/DISTlib/source/*.h src/mad_*.h)) +CC_SRC :=$(notdir $(wildcard libs/DISTlib/source/*.c src/mad_*.c)) + CC_HDR += madx.h SDDS.h CC_SRC += + + # files specific dependency flags $(OBJDIR)/mad_gxx11c.d: CPPFLAGS += $(call incdir,/usr/X11/include /opt/X11/include) From 0ff4c06f5b4fde644a9544af026dc90a5209b47b Mon Sep 17 00:00:00 2001 From: Tobias Date: Mon, 23 Mar 2020 18:57:58 +0100 Subject: [PATCH 04/20] Forgotten file. --- src/mad_dist.c | 140 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 140 insertions(+) create mode 100644 src/mad_dist.c diff --git a/src/mad_dist.c b/src/mad_dist.c new file mode 100644 index 000000000..a706642dc --- /dev/null +++ b/src/mad_dist.c @@ -0,0 +1,140 @@ +#include "madx.h" +void setdistparameters(char *type, int cut_l, double* cuts, int index, int *dist_type, double *start_v, double *stop_v, double *coord_t){ + + //1st First input (coordinate) + //2nd coord type, 0 - action angle, 1 - normalized + //3rd type 4 - uniform dist, 5 - normalized, 6-reyligh + //5 start value for uniform, mean value for normalized and reyligh + //6 sstop value for uniform, sigma for normalized and reyligh + + if(strcmp(type, "gauss")==0 ){ + dist_type[index] = 6; + dist_type[index+1] = 4; + start_v[index] = 0; + start_v[index+1] = 0; + stop_v[index] = 1; + stop_v[index+1] = twopi; + for(int i=0; i < 6; i++) + coord_t[i] = 0; + + if(cut_l==2){ + setactionanglecut(index, cuts[0], cuts[1]); + } + } + else if(strcmp(type, "uniform")==0){ + dist_type[index] = 4; + dist_type[index+1] = 4; + start_v[index] = pow(cuts[0],2); + start_v[index+1] = 0; + stop_v[index] = pow(cuts[1],2); + stop_v[index+1] = twopi; + for(int i=0; i < 6; i++) + coord_t[i] = 0; + } + else if(strcmp(type, "fixed")==0){ + dist_type[index] = 0; + dist_type[index+1] = 0; + start_v[index] = 2; //2J + start_v[index+1] = 0; + stop_v[index] = 0; + stop_v[index+1] = 0; + for(int i=0; i < 6; i++) + coord_t[i] = 0; + } + else if(strcmp(type, "zero")==0){ + dist_type[index] = 0; + dist_type[index+1] = 0; + start_v[index] = 0; //2J + start_v[index+1] = 0; + stop_v[index] = 0; + stop_v[index+1] = 0; + for(int i=0; i < 6; i++) + coord_t[i] = 0; + } +} + + +void pro_distribution(struct in_cmd* p){ + struct table* dist_t; + char *table_name; + char *type; + double num ; + double eigen [36]; + int dist_type [6]; + double start_v[6], stop_v[6], coord_t[6]; + int cut_l; + double cuts[3]; + int const C_HOR = 0,C_VER = 2,C_LON = 4; + + double npart_d = command_par_value("npart", p->clone); + int npart = (int)npart_d; + double x[npart], px[npart], y[npart], py[npart], t[npart], pt[npart]; + + initializedistribution(1); + setemitt12(get_value("beam", "ex"), get_value("beam", "ey")); + setemitt3(get_value("beam", "et")); + sete0andmass0(get_value("beam", "energy"),get_value("beam", "mass") ); + trbegn_(oneturnmat,eigen); + settasmatrixtranspose(eigen); + settotalsteps(npart); + + type = command_par_string("horizontal", p->clone); + cut_l = command_par_vector("cutsig_h", p->clone, cuts); + setdistparameters(type, cut_l, cuts, C_HOR, dist_type, start_v,stop_v, coord_t); + + type = command_par_string("vertical", p->clone); + cut_l = command_par_vector("cutsig_v", p->clone, cuts); + setdistparameters(type, cut_l, cuts, C_VER, dist_type, start_v,stop_v, coord_t); + + type = command_par_string("longitudinal", p->clone); + cut_l = command_par_vector("cutsig_l", p->clone, cuts); + setdistparameters(type, cut_l, cuts, C_LON, dist_type, start_v,stop_v, coord_t); + + + for(int i=0; i<6; i++){ + setscan_para_diagonal(i,coord_t[i],dist_type[i],start_v[i],stop_v[i]); + } + getunconvertedcoord(x,px,y,py,t,pt, &npart); + + + table_name = command_par_string("table", p->clone); + dist_t = make_table(table_name, "Distribution", dist_table_cols,dist_table_types, npart); + add_to_table_list(dist_t, table_register); + + + for(int i=0; i< npart; i++){ + double_to_table_curr(table_name, "x", &x[i]); + double_to_table_curr(table_name, "px", &px[i]); + double_to_table_curr(table_name, "y", &y[i]); + double_to_table_curr(table_name, "py", &py[i]); + double_to_table_curr(table_name, "t", &t[i]); + double_to_table_curr(table_name, "pt", &pt[i]); + num = (int)i; + double_to_table_curr(table_name, "number", &num); + augmentcountonly(table_name); + } + + if (dist_t->header == NULL) dist_t->header = new_char_p_array(40); + + //strncpy(tmp, t->org_sequ->name, NAME_L); + //table_add_header(t, "@ SEQUENCE %%%02ds \"%s\"", strlen(tmp),stoupper(tmp)); + char tmp[10]; + int i = get_string("beam", "particle", tmp); + + table_add_header(dist_t, "@ PARTICLE %%%02ds \"%s\"", i, stoupper(tmp)); + table_add_header(dist_t, "@ MASS %%le %F", get_value("beam", "mass")); + table_add_header(dist_t, "@ CHARGE %%le %F", get_value("beam", "charge")); + table_add_header(dist_t, "@ ENERGY %%le %F", get_value("beam", "energy")); + table_add_header(dist_t, "@ PC %%le %F", get_value("beam", "pc")); + table_add_header(dist_t, "@ EX %%le %F", get_value("beam", "ex")); + table_add_header(dist_t, "@ EY %%le %F", get_value("beam", "ey")); + table_add_header(dist_t, "@ ET %%le %F", get_value("beam", "et")); + type = command_par_string("horizontal", p->clone); + table_add_header(dist_t, "@ DIST_TYPE_H %%%02ds \"%s\"", strlen(type), stoupper(type)); + type = command_par_string("vertical", p->clone); + table_add_header(dist_t, "@ DIST_TYPE_V %%%02ds \"%s\"", strlen(type), stoupper(type)); + type = command_par_string("longitudinal", p->clone); + table_add_header(dist_t, "@ DIST_TYPE_L %%%02ds \"%s\"", strlen(type), stoupper(type)); + + free_distribution(0); // free the memory in the distribution block +} \ No newline at end of file From 8018b3c1421d88a6cb989b48b071223d73debcb9 Mon Sep 17 00:00:00 2001 From: Tobias Date: Mon, 23 Mar 2020 23:28:09 +0100 Subject: [PATCH 05/20] Another file. --- src/mad_dist.h | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 src/mad_dist.h diff --git a/src/mad_dist.h b/src/mad_dist.h new file mode 100644 index 000000000..3d4c65466 --- /dev/null +++ b/src/mad_dist.h @@ -0,0 +1,8 @@ +#ifndef MAD_DIST_H +#define MAD_DIST_H + + +void pro_distribution(struct in_cmd* p); +void setdistparameters(char *type, int cut_l, double* cuts, int index, int *dist_type, double *start_v, double *stop_v, double *coord_t); + +#endif // MAD_DIST_H \ No newline at end of file From dc8dcc8a5ee673a270131d07dc1bd0c2b8c0b7cb Mon Sep 17 00:00:00 2001 From: Tobias Date: Tue, 24 Mar 2020 12:13:43 +0100 Subject: [PATCH 06/20] Works with old compilation. --- libs/DISTlib/source/distgeneration.c | 46 ++- libs/DISTlib/source/distgeneration.h | 16 +- libs/DISTlib/source/distinterface.c | 19 +- libs/DISTlib/source/distinterface.h | 10 +- libs/DISTlib/source/file_reader.c | 526 --------------------------- libs/DISTlib/source/file_reader.h | 18 - libs/DISTlib/source/helper.c | 7 +- libs/DISTlib/source/helper.h | 7 +- 8 files changed, 70 insertions(+), 579 deletions(-) delete mode 100644 libs/DISTlib/source/file_reader.c delete mode 100644 libs/DISTlib/source/file_reader.h diff --git a/libs/DISTlib/source/distgeneration.c b/libs/DISTlib/source/distgeneration.c index fcc9abf08..b6936d46c 100644 --- a/libs/DISTlib/source/distgeneration.c +++ b/libs/DISTlib/source/distgeneration.c @@ -6,8 +6,6 @@ #include "distgeneration.h" #include "distinterface.h" #include "outputdist.h" -#include "file_reader.h" - void gensixcanonical(){ @@ -72,9 +70,6 @@ int gettotalgridlength(){ void generate_grid(){ int counter = 0; - double tc[6]; - double tmp[6]; - double tmp_n[6]; double **readin; int totallenght = gettotalgridlength(); readin = (double**)malloc(dim*sizeof(double*)); @@ -145,7 +140,6 @@ void normalized2canonical(double normalized_in[6], double cancoord[6]){ double deltap = normalized[4]; double deltas = normalized[5]; double *xap; - double det = (dist->tas[4][4]*dist->tas[5][5] - dist->tas[4][5]*dist->tas[5][4]); for(int i=0; i<4;i++){ lindeltas = lindeltas+dist->tas[4][i]*normalized[i]; lindp=lindp+dist->tas[5][i]*normalized[i]; @@ -289,8 +283,8 @@ double randn(double mu, double sigma) do { - U1 = -1 + ((double) rand () / RAND_MAX) * 2; - U2 = -1 + ((double) rand () / RAND_MAX) * 2; + U1 = -1 + ( rand () / (double)RAND_MAX) * 2; + U2 = -1 + ( rand () / (double)RAND_MAX) * 2; W = pow (U1, 2) + pow (U2, 2); } while (W >= 1 || W == 0); @@ -307,11 +301,45 @@ double randn(double mu, double sigma) double rand_uni(double low, double high) { - return ( (double)rand() * ( high - low ) ) / (double)RAND_MAX + low; + return (double)(rand() * ( high - low ) ) / (double)RAND_MAX + low; } double randray(double mu, double sigma){ double low = 0; double high =1; return pow((mu+sigma*sqrt((-2*log(rand_uni(low, high))))),2); +} + + +void allocateincoord(int linecount){ + dist->incoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); + dist->outcoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); + dist->totincoord = linecount; + for(int i=0; i< dim; i++){ + dist->ref->typeused[i]=-1; + } + for(int i=0; iincoord[i] = (struct coordinates*)malloc(sizeof(struct coordinates)); + + dist->incoord[i]->coord = (double*)malloc(dim*sizeof(double)); + dist->incoord[i]->readin = (double*)malloc(32*sizeof(double)); + + dist->incoord[i]->mass=0; + dist->incoord[i]->a=0; + dist->incoord[i]->z=0; + + dist->outcoord[i] = (struct coordinates*)malloc(sizeof(struct coordinates)); + dist->outcoord[i]->coord = (double*)malloc(dim*sizeof(double)); + + } + dist->isallocated =1; +} +void deallocateincoord(void){ + + + free(dist->incoord); + free(dist->outcoord); + //dist->ref->typeused = (int*)malloc(dim*sizeof(int)); + dist->isallocated =0; + dist->totincoord = -1; } \ No newline at end of file diff --git a/libs/DISTlib/source/distgeneration.h b/libs/DISTlib/source/distgeneration.h index df6d6aa22..2d2da4939 100644 --- a/libs/DISTlib/source/distgeneration.h +++ b/libs/DISTlib/source/distgeneration.h @@ -76,14 +76,14 @@ struct appliedcut{ //void action2sixinternal_(double tc[6], double *results, double *normalized); //int checkdist(); //void createrandomdist_(); -void gensixcanonical(); +void gensixcanonical(void); int particle_within_limits_physical(double *physical); int particle_within_limits_normalized(double *normalized); int particle_with_limits_action(int i, double value); -void generatefromnormalized(); -void generatefrommixed(); -void generatefromaction(); -void generatefromphysical(); +void generatefromnormalized(void); +void generatefrommixed(void); +void generatefromaction(void); +void generatefromphysical(void); void action2normalized(double acangl[6], double normalized[6]); void normalized2canonical(double normalized[6], double cancoord[6]); double randn(double mu, double sigma); @@ -91,5 +91,7 @@ double rand_uni(double low, double high); void createLinearSpaced(int length, double start, double stop, double *eqspaced ); double randray(double mu, double sigma); void createcoordinates(int index, double start, double stop, int length, int type); -int gettotalgridlength(); -void generate_grid(); \ No newline at end of file +int gettotalgridlength(void); +void generate_grid(void); +void allocateincoord(int linecount); +void deallocateincoord(void); \ No newline at end of file diff --git a/libs/DISTlib/source/distinterface.c b/libs/DISTlib/source/distinterface.c index bb277618a..95a25f1ff 100644 --- a/libs/DISTlib/source/distinterface.c +++ b/libs/DISTlib/source/distinterface.c @@ -6,7 +6,6 @@ #include "distinterface.h" #include "distgeneration.h" #include "outputdist.h" -#include "file_reader.h" /* @@ -241,7 +240,7 @@ void get6trackcoord(double *x, double *xp, double *y, double *yp, double *sigma, nparticles = *totparticles; for(int i=0; i < nparticles; i++){ - canonical2six(dist->outcoord[i]->coord, dist->ref->beta0, dist->ref->pc0, dist->ref->mass0, dist->incoord[i]->mass, tmp); + canonical2six(dist->outcoord[i]->coord, dist->ref->beta0, dist->ref->pc0, dist->incoord[i]->mass, tmp); x[i] = tmp[0]; xp[i] = tmp[1]; y[i] = tmp[2]; @@ -253,7 +252,6 @@ void get6trackcoord(double *x, double *xp, double *y, double *yp, double *sigma, } void getunconvertedcoord(double *x, double *xp, double *y, double *yp, double *sigma, double *deltap, int *totparticles){ - double tmp[6]; int nparticles; if(dist->isDistrcalculated ==0){ gensixcanonical(); @@ -296,19 +294,20 @@ void readtasmatrixfile(const char* filename_in){ void getrefpara(double *energy0, double *mass0, int *a0, int *z0){ - *energy0=dist->ref->e0; - *mass0=dist->ref->mass0; - *a0=dist->ref->a0; - *z0=dist->ref->z0; + *energy0= dist->ref->e0; + *mass0 = dist->ref->mass0; + *a0 = dist->ref->a0; + *z0 = dist->ref->z0; } -int readfile_f(const char* filename_in, int strlen){ +/* +void readfile_f(const char* filename_in, int strlen){ char filename [strlen]; strncpy(filename, filename_in, strlen); filename[strlen] = '\0'; readfile(filename); -} +}*/ -int writefile_f(const char* filename_in, int strlen){ +void writefile_f(const char* filename_in, int strlen){ char filename [strlen]; strncpy(filename, filename_in, strlen); filename[strlen] = '\0'; diff --git a/libs/DISTlib/source/distinterface.h b/libs/DISTlib/source/distinterface.h index d2a9756ca..f0e05a676 100644 --- a/libs/DISTlib/source/distinterface.h +++ b/libs/DISTlib/source/distinterface.h @@ -9,8 +9,8 @@ struct distparam* dist; struct distparam* diststart; int dim; int distn; -int writefile_f(const char* filename_in, int strlen); -int readfile_f(const char* filename_in, int strlen); +void writefile_f(const char* filename_in, int strlen); +void readfile_f(const char* filename_in, int strlen); void readtasmatrixfile(const char* filename_in); void initializedistribution(int numberOfDist); void setdistribution(int ndist); @@ -31,4 +31,8 @@ void setdisptas(double dx, double dpx, double dy, double dpy); void settwisstas(double betax, double alfax, double betay, double alfay); void getunconvertedcoord(double *x, double *xp, double *y, double *yp, double *sigma, double *deltap, int *totparticles); void setactionanglecut(int variable, double min, double max); -void free_distribution(int i); \ No newline at end of file +void free_distribution(int i); +void settasmatrix_element(double value, int row, int column); +void addclosedorbit(double *clo); +void getarraylength(int *totlength); +void getrefpara(double *energy0, double *mass0, int *a0, int *z0); \ No newline at end of file diff --git a/libs/DISTlib/source/file_reader.c b/libs/DISTlib/source/file_reader.c deleted file mode 100644 index 21966e6e7..000000000 --- a/libs/DISTlib/source/file_reader.c +++ /dev/null @@ -1,526 +0,0 @@ -#include -#include -#include -#include -#include -#include "helper.h" -#include "distinterface.h" -#include "outputdist.h" -#include "file_reader.h" -#include "distgeneration.h" -//#include "round_near.c" - - -int readfile(const char* filename){ - dist->incoordtype=2; - FILE *file = fopen ( filename, "r" ); - double** table = malloc(MAX_ROWS * sizeof(double*)); // allocate the rows - for (int r = 0; r < MAX_ROWS; ++r){ - table[r] = malloc(MAX_COLUMNS * sizeof(double)); // allocate the columns - } - - int linecount = -1; - int numcolum, index; - char columns[MAX_COLUMNS][MAX_LENGTH]; - char units[MAX_COLUMNS][MAX_LENGTH]; - char line [ MAX_LENGTH*50 ]; /* or other suitable maximum line size */ - char tosplit [ MAX_LENGTH*10 ]; - char value_s[MAX_LENGTH], shorty[MAX_LENGTH]; - char * parameter = (char*)malloc(MAX_LENGTH*sizeof(char)); - char * parameter_tmp = (char*)malloc(MAX_LENGTH*sizeof(char)); - char * unit =(char*)malloc(MAX_LENGTH*sizeof(char)); - char * unit_tmp =(char*)malloc(MAX_LENGTH*sizeof(char)); - char * at =(char*)malloc(MAX_LENGTH*sizeof(char)); - double value; - char *e ; - if ( file != NULL ){ - while ( fgets ( line, sizeof line, file ) != NULL ){ /* read a line */ - double multifactor = 1; - int k =0; - if(strncmp(line, "@", 1)==0){ - while( line[k] ) { - line[k] = (tolower(line[k])); - k++; - } - sscanf( line, "%s %s %s", at, parameter_tmp, value_s); - if(strncmp(value_s, "%", 1)!=0){ - unit_tmp =strstr(parameter_tmp, "["); - - if(unit_tmp != NULL){ - strcpy(unit,unit_tmp); - e = strchr(parameter_tmp, '['); - index = (int)(e - parameter_tmp); - strncpy(parameter, parameter_tmp, index); - } - else{ - strcpy(unit, "nounit"); - strcpy(parameter, parameter_tmp); - } - - value = strtod(value_s,NULL); - - if(strcmp(parameter, "energy0")==0){ - multifactor = getEnergyUnit(unit); - dist->ref->e0=value*multifactor; - } - else if(strcmp(parameter, "pc0")==0){ - multifactor = getEnergyUnit(unit); - dist->ref->pc0=value*multifactor; - } - else if(strcmp(parameter, "a0")==0){ - dist->ref->a0=round(value); - } - else if(strcmp(parameter, "z0")==0){ - dist->ref->z0=round(value); - } - else if(strcmp(parameter, "mass0")==0){ - multifactor = getEnergyUnit(unit); - dist->ref->mass0=value*multifactor; - } - else if(strcmp(parameter, "charge0")==0){ - dist->ref->charge0=round(value); - } - else{ - issue_error2("Not recogniced parameter!", parameter); - } - } - } - else if(strncmp(line, "!", 1)==0 || strncmp(line, "$", 1)==0){ - printf("Comment line %s \n", line ); - } - else if(strncmp(line, "!", 1)!=0 && linecount==-1){ - - strcpy(tosplit, line); - numcolum = splitline(tosplit, columns, units); - linecount++; - } - else if(linecount>=0 && strncmp(line, "$", 1)!=0){ - add2table(table, line, linecount); - linecount++; - } - else{ - issue_error2("Something went wrong with reading line", line); - } - - } - fclose ( file ); - } - else - { - perror ( filename ); - } - allocateincoord(linecount); //alocates the memory - fillcoordstructure(numcolum,columns, units, table); //fils the coordinate structures - calculaterefparam(); //Sets up the necessary ref parameters and checks for consistensy. - convert2standard(); // converts to the units used internal , x, px, y, py, zeta, deltap - checkinputtype(); - free_readin_memory(table); - return 0; -} -void checkinputtype(){ - if(dist->ref->typeused[0]==dist->ref->typeused[1]==dist->ref->typeused[2]==dist->ref->typeused[3]==dist->ref->typeused[4]==dist->ref->typeused[5]) - dist->incoordtype=dist->ref->typeused[0]; - else if(dist->ref->typeused[4]==2==dist->ref->typeused[5]==2 && dist->ref->typeused[0]==dist->ref->typeused[1]==dist->ref->typeused[2]==dist->ref->typeused[3]==0) - dist->incoordtype=3; //Mixed with action angle coordinates - else if(dist->ref->typeused[4]==2==dist->ref->typeused[5]==2 && dist->ref->typeused[0]==dist->ref->typeused[1]==dist->ref->typeused[2]==dist->ref->typeused[3]==1) - dist->incoordtype=4; //Mixed with normalised coordinates - else - issue_error("Non comptable types have been entered."); - - if(dist->incoordtype==0){ - - } - else if(dist->incoordtype==1){ - - } - else if(dist->incoordtype==3){ - - } - else if(dist->incoordtype==4){ - - } - -} - -// This checks that reference energy is only set once and sets up the mass in case it is not a column -void calculaterefparam(){ - if(dist->ref->mass0 == 0){ - issue_error("A mass is needed! \n"); - } - if(dist->ref->pc0 == 0 && dist->ref->e0==0){ - issue_error("A energy is needed! \n"); - } - if(dist->ref->pc0 > 0 && dist->ref->e0 > 0){ - issue_error("Can't set both energy and momentum! \n"); - } - if(dist->ref->pc0 > 0){ - dist->ref->e0 = momentum2energy(dist->ref->pc0,dist->ref->mass0); - } - if(dist->ref->e0 > 0){ - dist->ref->pc0 = energy2momentum(dist->ref->e0,dist->ref->mass0); - } - - dist->ref->beta0 = (dist->ref->pc0)/(dist->ref->e0); - - if(dist->incoord[0]->mass < 1e-16){ //enough to check the first - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->mass = dist->ref->mass0; //Gives all the same mass. - } - } - if(dist->incoord[0]->a ==0){ //enough to check the first - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->a = dist->ref->a0; //Gives all the same mass. - } - } - if(abs(dist->incoord[0]->z) ==0){ //enough to check the first - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->z = dist->ref->z0; //Gives all the same mass. - } - } -} - -//Converts to the internal used canonical -void convert2standard(){ - for(int i=0; itotincoord; i++){ - for(int j=0; j<6; j++){ - dist->incoord[i]->coord[j]= dist->incoord[i]->readin[j]; // have to have checks here..! - } - } - - if(dist->ref->en_like==-1){ - issue_info("No energy variable is set. Assume 0 deviation from reference energy \n"); - } - else if(dist->ref->en_like==0){ //energy - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->coord[5] = (momentum2energy(dist->incoord[i]->readin[6], dist->incoord[i]->mass)-(dist->ref->pc0))/(dist->ref->pc0); - } - } - else if(dist->ref->en_like==1){ //momentum - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->coord[5] = ((dist->incoord[i]->readin[7], dist->incoord[i]->mass)-(dist->ref->pc0))/(dist->ref->pc0); - } - } - else if(dist->ref->en_like==2){ //psigma - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->coord[5] = psigma2deltap(dist->incoord[i]->readin[8], dist->ref->beta0); - } - } - else if(dist->ref->en_like==3){ //pt - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->coord[5] = pt2deltap(dist->incoord[i]->readin[9], dist->ref->beta0); - } - } - - if(dist->ref->time_like==-1){ - issue_info("No time variable is set. Assume 0 deviation \n"); - } - else if(dist->ref->time_like==2){ //sigma - double beta, pc; - for(int i=0; i< dist->totincoord; i++){ - pc = dist->ref->pc0+dist->incoord[i]->coord[5]; - beta = (pc)/momentum2energy(pc, dist->incoord[i]->mass); - dist->incoord[i]->coord[4] = sigma2zeta(dist->incoord[i]->readin[10], dist->ref->beta0, beta); - } - } - else if(dist->ref->time_like==3){ //t or tau - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->coord[4] = tau2zeta(dist->incoord[i]->readin[11], dist->ref->beta0); - } - } - if(dist->ref->ang_like==-1){ - issue_error("No angle variable is set. \n"); - } - else if(dist->ref->ang_like==1){ //t or tau - for(int i=0; i< dist->totincoord; i++){ - dist->incoord[i]->coord[1] = dist->incoord[i]->coord[6]/(1+dist->incoord[i]->coord[5]); - dist->incoord[i]->coord[3] = dist->incoord[i]->coord[7]/(1+dist->incoord[i]->coord[5]); - } - } -} - -void allocateincoord(int linecount){ - dist->incoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); - dist->outcoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); - dist->totincoord = linecount; - for(int i=0; i< dim; i++){ - dist->ref->typeused[i]=-1; - } - for(int i=0; iincoord[i] = (struct coordinates*)malloc(sizeof(struct coordinates)); - - dist->incoord[i]->coord = (double*)malloc(dim*sizeof(double)); - dist->incoord[i]->readin = (double*)malloc(32*sizeof(double)); - - dist->incoord[i]->mass=0; - dist->incoord[i]->a=0; - dist->incoord[i]->z=0; - - dist->outcoord[i] = (struct coordinates*)malloc(sizeof(struct coordinates)); - dist->outcoord[i]->coord = (double*)malloc(dim*sizeof(double)); - - } - dist->isallocated =1; -} -void deallocateincoord(){ - - - for(int i=0; itotincoord; i++){ - - // free(dist->incoord[i]->coord); - // free(dist->incoord[i]->readin); - // free(dist->incoord[i]); - // free(dist->outcoord[i]->coord); -// free(dist->outcoord[i]); - - } - - free(dist->incoord); - free(dist->outcoord); - //dist->ref->typeused = (int*)malloc(dim*sizeof(int)); - dist->isallocated =0; - dist->totincoord = -1; -} - - - -void setreadin(int coordorder, int column, double ** table, double multifactor, int realorder){ - if(dist->ref->typeused[realorder]!=-1){ - issue_error("Only allowed to set one coordinates index ones"); - } - - if(coordorder < 20){ - dist->ref->typeused[realorder] = 2; - } - else if(coordorder >= 20 <=25){ - dist->ref->typeused[realorder] = 1; - } - else{ - dist->ref->typeused[realorder] = 0; - } - for(int i=0;i < dist->totincoord; i++){ - dist->incoord[i]->readin[coordorder] = multifactor*table[i][column]; - } -} - - -int splitline(char* line, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH] ){ - char* word; - int count = 0; - int k=0; - - while( line[k] ) { - line[k] = (tolower(line[k])); - k++; - } - - word = strtok(line, " "); - if(strncmp("*", word, 1)!=0){ - add2internaltab(word, columns, units, count); - count++; - } - while ((word = strtok(NULL, " ")) != NULL){ - add2internaltab(word, columns, units, count); - count ++; - } - - return count; -} -void add2internaltab(char* word, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH] , int count){ - - - char * unit =(char*)malloc(MAX_LENGTH*sizeof(char)); - char * word_tmp =(char*)malloc(MAX_LENGTH*sizeof(char)); - char * columnname =(char*)malloc(MAX_LENGTH*sizeof(char)); - char *e; - int index; - strcpy(word_tmp, word); - unit =strstr(word_tmp, "["); - if(unit != NULL){ - - strcpy(units[count], unit); - e = strchr(word_tmp, '['); - index = (int)(e - word_tmp); - strncpy(columns[count], word_tmp, index); - columns[count][index]='\0'; - } - else{ - strcpy(units[count], "nounit"); - strcpy(columns[count], word_tmp); - } -} - -void add2table(double ** table, char* line, int linenum){ - char* word; - int count = 0; - word = strtok(line, " "); - table[linenum][0] = atof(word); - count++; - - /* the following loop gets the rest of the words until the - * end of the message */ - while ((word = strtok(NULL, " ")) != NULL){ - table[linenum][count] = atof(word); - count ++; - } -} - -void fillcoordstructure(int numcolum, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH], double **table){ - - double multifactor; - for(int i=0; i< numcolum; i++){ - if(strcmp(columns[i], "x")==0){ - multifactor = getMetricUnit(units[i]); - setreadin(0, i, table, multifactor, 0); - } - else if(strcmp(columns[i], "px")==0){ - checkifangset(0); - setreadin(1, i, table,1, 1); - } - else if(strcmp(columns[i], "y")==0){ - multifactor = getMetricUnit(units[i]); - setreadin(2, i, table, multifactor, 2); - } - else if(strcmp(columns[i], "py")==0){ - checkifangset(0); - setreadin(3, i, table,1, 3); - } - else if(strcmp(columns[i], "zeta")==0){ - checkiftimeset(5); - multifactor = getMetricUnit(units[i]); - setreadin(4, i, table, multifactor, 4); - } - else if(strcmpnl(columns[i], "deltap")==0){ - checkifenergyset(5); - setreadin(5, i, table, 1, 5); - } - else if(strcmp(columns[i], "energy")==0){ - multifactor = getEnergyUnit(units[i]); - checkifenergyset(0); - //setnonstandard(0, i, table); - setreadin(6, i, table, 1, 5); - } - else if(strcmp(columns[i], "pc")==0){ - multifactor = getEnergyUnit(units[i]); - checkifenergyset(1); - //setnonstandard(1, i, table); - setreadin(7, i, table, 1, 5); - } - else if(strcmp(columns[i], "psigma")==0){ - checkifenergyset(2); - // setnonstandard(2, i, table); - setreadin(8, i, table, 1, 5); - } - else if(strcmp(columns[i], "ptau")==0 || strcmp(columns[i], "pt")==0 ){ - checkifenergyset(3); - // setnonstandard(3, i, table); - setreadin(9, i, table, 1, 4); - } - else if(strcmp(columns[i], "sigma")==0){ - checkiftimeset(2); - multifactor = getMetricUnit(units[i]); - //setnonstandard(4, i, table); - setreadin(10, i, table, 1, 4); - } - else if(strcmp(columns[i], "tau")==0 || strcmp(columns[i], "t")==0 ){ - checkiftimeset(3); - multifactor = getMetricUnit(units[i]); - //setnonstandard(5, i, table); - setreadin(11, i, table, 1, 4); - } - else if(strcmp(columns[i], "xp")==0){ - checkifangset(1); - //setnonstandard(6, i, table); - setreadin(12, i, table, 1, 1); - } - else if(strcmp(columns[i], "yp")==0){ - checkifangset(1); - //setnonstandard(7, i, table); - setreadin(13, i, table, 1, 3); - } - else if(strcmp(columns[i], "jx")==0) - setreadin(20, i, table, 1, 0); - else if(strcmp(columns[i], "phix")==0) - setreadin(21, i, table, 1, 1); - else if(strcmp(columns[i], "jy")==0) - setreadin(22, i, table, 1, 2); - else if(strcmp(columns[i], "phiy")==0) - setreadin(23, i, table, 1, 3); - else if(strcmp(columns[i], "jz")==0) - setreadin(24, i, table, 1, 4); - else if(strcmp(columns[i], "phiz")==0) - setreadin(25, i, table, 1, 5); - else if(strcmp(columns[i], "xn")==0) - setreadin(26, i, table, 1, 0); - else if(strcmp(columns[i], "pxn")==0) - setreadin(27, i, table, 1, 1); - else if(strcmp(columns[i], "yn")==0) - setreadin(28, i, table, 1, 2); - else if(strcmp(columns[i], "pyn")==0) - setreadin(29, i, table, 1, 3); - else if(strcmp(columns[i], "zn")==0) - setreadin(30, i, table, 1, 4); - else if(strcmp(columns[i], "pzn")==0) - setreadin(31, i, table, 1, 5); - else if(strcmp(columns[i], "mass")==0){ - multifactor = getEnergyUnit(units[i]); - for(int j=0;j < dist->totincoord; j++){ - dist->incoord[j]->mass = multifactor*table[j][i]; - } - } - else if(strcmp(columns[i], "a")==0){ - for(int j=0;j < dist->totincoord; j++){ - dist->incoord[j]->a = table[j][i]; - } - } - else if(strcmp(columns[i], "z")==0){ - for(int j=0;j < dist->totincoord; j++){ - dist->incoord[j]->z = table[j][i]; - } - } - } -} -double getEnergyUnit(char *unit){ - double multif; - if(strcmp(unit, "[ev]")==0){ - multif = 1e-9; - } - else if(strcmp(unit, "[kev]")==0){ - multif = 1e-6; - } - else if(strcmp(unit, "[mev]")==0){ - multif = 1e-3; - } - else if(strcmp(unit, "[gev]")==0 || strcmp(unit, "nounit")==0){ - multif = 1; - } - else if(strcmp(unit, "[tev]")==0){ - multif = 1e3; - } - else{ - issue_error("Unknow unit. Must be of the form [Mev]"); - } - return multif; -} - -double getMetricUnit(char *unit){ - double multif; - - if(strcmp(unit, "[mm]")==0){ - multif = 1e-3; - } - else if(strcmp(unit, "[m]")==0 || strcmp(unit, "nounit")==0 ){ - multif = 1; - } - else{ - issue_warning("Unknow unit. Must be of the form [Mev]"); - } - return multif; -} - -void free_readin_memory(double **table){ - - for (int r = 0; r < MAX_ROWS; ++r){ - free(table[r]); - } - free(table); -} \ No newline at end of file diff --git a/libs/DISTlib/source/file_reader.h b/libs/DISTlib/source/file_reader.h deleted file mode 100644 index de5cda18c..000000000 --- a/libs/DISTlib/source/file_reader.h +++ /dev/null @@ -1,18 +0,0 @@ -int readfile(const char* filenamet); -int splitline(char* line, char columns[MAX_COLUMNS][MAX_LENGTH],char units[MAX_COLUMNS][MAX_LENGTH] ); -void add2table(double ** table, char* line, int linenum ); -void allocateincoord(int linecount); -void setreadin(int coordorder, int column, double ** table, double multifactor, int realcoord); -void checkifenergyset(int entype); -void checkiftimeset(int entype); -void checkifangset(int entype); -void calculaterefparam(); -void convert2standard(); -void add2internaltab(char* word, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH], int count) ; -double getEnergyUnit(char *unit); -double getMetricUnit(char *unit); -void fillcoordstructure(int numcolum, char columns[MAX_COLUMNS][MAX_LENGTH], char units[MAX_COLUMNS][MAX_LENGTH], double **table); -void free_readin_memory(double **table); -void checkinputtype(); -void deallocateincoord(); -void allocateincoord(); \ No newline at end of file diff --git a/libs/DISTlib/source/helper.c b/libs/DISTlib/source/helper.c index 6b28b653e..7504f0954 100644 --- a/libs/DISTlib/source/helper.c +++ b/libs/DISTlib/source/helper.c @@ -9,7 +9,7 @@ #include "outputdist.h" /*This function converts from canoncial to sixtrack tracking variables*/ -void canonical2six(double *canonical, double beta0, double pc0, double mass0, double mass, double *coord){ +void canonical2six(double *canonical, double beta0, double pc0, double mass, double *coord){ double deltap = canonical[5]; double beta = (pc0+deltap)/momentum2energy((pc0+deltap), mass); double rv = beta0/beta; @@ -126,8 +126,7 @@ void solve2by2eq(double a1, double b1, double c1, double a2, double b2, double c void mtrx_vector_mult_pointer(int mp, int np, double **mtrx_a, double mtrx_b[6], double result[6]) { - int m = mp; - int n = np; + result[0]=0; result[1]=0; result[2]=0; @@ -135,7 +134,7 @@ void mtrx_vector_mult_pointer(int mp, int np, double **mtrx_a, double mtrx_b[6] result[4]=0; result[5]=0; - register int i=0, j=0, k=0; + register int i=0, k=0; for (i = 0; i < mp; i++) { for (k = 0; k < np; k++) diff --git a/libs/DISTlib/source/helper.h b/libs/DISTlib/source/helper.h index bb189d4b2..026b6f99f 100644 --- a/libs/DISTlib/source/helper.h +++ b/libs/DISTlib/source/helper.h @@ -1,7 +1,7 @@ #define MAX_COLUMNS 100 #define MAX_LENGTH 100 #define MAX_ROWS 100000 -void canonical2six(double *canonical, double beta0, double pc0, double mass0, double mass, double *coord); +void canonical2six(double *canonical, double beta0, double pc0, double mass, double *coord); double momentum2energy(double momentum, double mass); double energy2momentum(double energy, double mass); double pt2deltap(double pt, double beta0 ); @@ -14,4 +14,7 @@ int strcmpnl (const char *s1, const char *s2); double sigma2zeta(double sigma, double beta0, double beta); double tau2zeta(double tau, double beta); void solve2by2eq(double a1, double b1, double c1, double a2, double b2, double c2, double *x); -void mtrx_vector_mult_pointer(int mp, int np, double **mtrx_a, double mtrx_b[6], double result[6]); \ No newline at end of file +void mtrx_vector_mult_pointer(int mp, int np, double **mtrx_a, double mtrx_b[6], double result[6]); +void checkiftimeset(int entype); +void checkifangset(int entype); +void checkifenergyset(int entype); \ No newline at end of file From f12bd75b1dc6ae371dce166067a981140de2e35c Mon Sep 17 00:00:00 2001 From: Tobias Date: Tue, 24 Mar 2020 13:35:18 +0100 Subject: [PATCH 07/20] Remove print. --- libs/DISTlib/source/distgeneration.c | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/libs/DISTlib/source/distgeneration.c b/libs/DISTlib/source/distgeneration.c index b6936d46c..544e85c4d 100644 --- a/libs/DISTlib/source/distgeneration.c +++ b/libs/DISTlib/source/distgeneration.c @@ -176,7 +176,6 @@ int particle_within_limits_physical(double *physical){ int particle_with_limits_action(int i, double value){ if(dist->cuts2apply->action[i]->isset==1){ - printf("iiiiiii %d \n", dist->cuts2apply->action[i]->isset); if(value > pow(dist->cuts2apply->action[i]->min,2) && value < pow(dist->cuts2apply->action[i]->max,2)) return 1; else return 0; } @@ -205,7 +204,6 @@ void createcoordinates(int index, double start, double stop, int length, int ty for(int i=0;i incoord[i]->coord[index] = start; } - printf("tttttttttttttttt %d %f \n", index, start); return; } @@ -234,7 +232,6 @@ void createcoordinates(int index, double start, double stop, int length, int ty else if(type==4){ // uniform random for(int i=0;i incoord[i]->coord[index] = rand_uni(start, stop); - //printf("%f \n", dist->coord[index-1]->values[i] ); } } @@ -253,7 +250,7 @@ void createcoordinates(int index, double start, double stop, int length, int ty tmp = randray(start, stop); dist->incoord[i]->coord[index] = tmp; - printf("%f \n", sqrt(tmp/2)); + //printf("%f \n", sqrt(tmp/2)); } } else From 3de98a3b191e9d0c9a17126824a725dd3644529d Mon Sep 17 00:00:00 2001 From: Tobias Date: Tue, 24 Mar 2020 13:37:33 +0100 Subject: [PATCH 08/20] Remove out. --- src/mad_cmd.c | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mad_cmd.c b/src/mad_cmd.c index 6623a7c91..f28614085 100644 --- a/src/mad_cmd.c +++ b/src/mad_cmd.c @@ -180,8 +180,7 @@ exec_command(void) scan_in_cmd(p); /* match input command with clone + fill */ current_command = p->clone; - // main commands - printf("outtt %s \n", p->cmd_def->module ); + // main command_par_string_user if (strcmp(p->cmd_def->module, "control") == 0) control(p); else if (strcmp(p->cmd_def->module, "c6t") == 0) conv_sixtrack(p); else if (strcmp(p->cmd_def->module, "edit") == 0) seq_edit_main(p); From 71846cf0f8ce7386a447e415993fe83e7151624a Mon Sep 17 00:00:00 2001 From: Tobias Date: Wed, 13 May 2020 13:47:19 +0200 Subject: [PATCH 09/20] New command to change p0. --- src/mad_dict.c | 31 ++++++++++++++++++++ src/twiss.f90 | 77 ++++++++++++++++++++++++++++++++++++++++++++++---- src/util.f90 | 1 + 3 files changed, 104 insertions(+), 5 deletions(-) diff --git a/src/mad_dict.c b/src/mad_dict.c index ba013f1b5..56b5f7e45 100644 --- a/src/mad_dict.c +++ b/src/mad_dict.c @@ -3087,4 +3087,35 @@ const char *const_element_def = "slice = [i, 1], " "comments = [s, none, none]; " " " +"changefp0: element none 0 45 " +"at = [r, 1.e20], " +"l = [r, 0], " +"kmax = [r, 0], " +"kmin = [r, 0], " +"calib = [r, 0], " +"polarity = [r, 0], " +"tilt = [r, 0], " +"lrad = [r, 0], " +"magnet = [i, 0], " +"p0 = [r, -1], " +"apertype = [s, circle, circle], " +"aperture = [r, {0}], " +"aper_offset = [r, {0}], " +"aper_tol = [r, {0, 0, 0}], " +"aper_vx = [r, {-1}], " +"aper_vy = [r, {-1}], " +"slot_id = [i, none], " +"assembly_id = [i, none], " +"mech_sep = [r, 0], " +"v_pos = [r, 0], " +"model = [i, -1], " +"method = [i, -1], " +"exact = [i, -1, 1], " +"nst = [i, -1], " +"from = [s, none], " +"type = [s, none, none], " +"aper_tilt = [r, 0], " +"slice = [i, 1], " +"comments = [s, none, none]; " +" " ; diff --git a/src/twiss.f90 b/src/twiss.f90 index e59895b77..25b8889d0 100644 --- a/src/twiss.f90 +++ b/src/twiss.f90 @@ -3439,9 +3439,11 @@ SUBROUTINE tmmap(code,fsec,ftrk,orbit,fmap,ek,re,te,fcentre,dl) call tmbb(fsec,ftrk,orbit,fmap,re,te) case (code_marker) + ! nothing on purpose! case (code_wire) + call tmchenergy(ftrk,orbit,fmap,ek,re, te) ! nothing for now... case (code_dipedge) @@ -3464,7 +3466,8 @@ SUBROUTINE tmmap(code,fsec,ftrk,orbit,fmap,ek,re,te,fcentre,dl) case (code_rfmultipole) call tmrfmult(fsec,ftrk,orbit,fmap,ek,re,te) - + case (code_changerefpc) + call tmchenergy(ftrk,orbit,fmap,ek,re, te) case default !--- anything else: ! nil (23, 28, 34) @@ -4624,10 +4627,10 @@ SUBROUTINE tmmult(fsec,ftrk,orbit,fmap,re,te) di = (dr * y + di * x) / (iord) + f_errors(2*iord+1) dr = drt enddo - re(2,1) = - dr + ! re(2,1) = - dr re(2,3) = + di re(4,1) = + di - re(4,3) = + dr + ! re(4,3) = + dr endif !---- Add the missing focussing component of thin dipoles @@ -4648,6 +4651,7 @@ SUBROUTINE tmmult(fsec,ftrk,orbit,fmap,re,te) re(5,3) = - re(4,6) !---- Second-order terms (use X,Y from orbit tracking). + if (fsec) then if (nord .ge. 2) then dr = zero @@ -4669,7 +4673,7 @@ SUBROUTINE tmmult(fsec,ftrk,orbit,fmap,re,te) te(4,3,3) = - di endif endif - +te = zero end SUBROUTINE tmmult SUBROUTINE tmoct(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te) @@ -6318,6 +6322,7 @@ SUBROUTINE tmdrf(fsec,ftrk,orbit,fmap,dl,ek,re,te) fmap = dl .ne. zero !---- First-order terms. + re(1,2) = dl re(3,4) = dl re(5,6) = dl/(beta*gamma)**2 @@ -6334,12 +6339,66 @@ SUBROUTINE tmdrf(fsec,ftrk,orbit,fmap,dl,ek,re,te) te(5,4,4) = te(5,2,2) te(5,6,6) = te(1,2,6) * three / (beta * gamma) ** 2 endif - + !---- Track orbit. if (ftrk) call tmtrak(ek,re,te,orbit,orbit) end SUBROUTINE tmdrf +SUBROUTINE tmchenergy(ftrk,orbit,fmap,ek,re, te) + use twisslfi + use twiss_elpfi + use twissbeamfi, only : deltap, pc, gamma, energy, beta + use matrices, only : EYE + use math_constfi, only : zero, one + use phys_constfi, only : clight + implicit none + !----------------------------------------------------------------------* + ! Purpose: * + ! TRANSPORT map for change of reference Energy * + ! Input: * + ! ftrk (logical) if true, track orbit. * + ! Input/output: * + ! orbit(6) (double) closed orbit. * + ! Output: * + ! fmap (logical) if true, element has a map. * + ! ek(6) (double) kick due to element. * + ! re(6,6) (double) transfer matrix. * + !----------------------------------------------------------------------* + logical :: fsec, ftrk, fmap, fcentre, fringe + double precision :: orbit(6), ek(6), re(6,6), te(6,6,6) + + + double precision :: energy1, pc1, gamma1, pt, mass, beta1i + double precision, external :: node_value, get_value + integer, external :: el_par_vector + + mass = get_value('probe ','mass ') + re = EYE + ek = zero + te = zero + pt = orbit(6) + + !---- Transfer map. + fmap = .true. + + energy1 = pt*pc+energy + pc1 = sqrt(energy1**2-mass**2) + gamma1 = energy1/mass + beta1i = energy1/pc1 + + re(2,2) = pc/pc1 + re(4,4) = pc/pc1 + re(6,6) = pc/pc1 + + ek(6) = beta1i*(gamma/gamma1-one) + + if(ftrk) call tmtrak(ek,re,te,orbit,orbit) + + +end SUBROUTINE tmchenergy + + SUBROUTINE tmrf(fsec,ftrk,fcentre,orbit,fmap,el,ds,ek,re,te) use twisslfi use twiss_elpfi @@ -6370,6 +6429,7 @@ SUBROUTINE tmrf(fsec,ftrk,fcentre,orbit,fmap,el,ds,ek,re,te) double precision :: orbit(6), ek(6), re(6,6), te(6,6,6) double precision :: ek_f(6), re_f(6,6), te_f(6,6,6) double precision :: ek_tmp(6), re_tmp(6,6), te_tmp(6,6,6) + double precision :: ek_ch(6), re_ch(6,6), te_ch(6,6,6) integer :: elpar_vl double precision :: rfv, rff, rfl, dl, omega, vrf, phirf, bvk @@ -6406,6 +6466,7 @@ SUBROUTINE tmrf(fsec,ftrk,fcentre,orbit,fmap,el,ds,ek,re,te) rff = g_elpar(r_freq) rfl = g_elpar(r_lag) bvk = node_value('other_bv ') + !-- LD: 20.6.2014 (bvk=-1: not -V -> V but lag -> pi-lag) if (bvk .eq. -one) then @@ -6451,6 +6512,12 @@ SUBROUTINE tmrf(fsec,ftrk,fcentre,orbit,fmap,el,ds,ek,re,te) endif call tmcat(fsec,re,te,rw,tw,re,te) + + if (changeref_p0) then + call tmchenergy(ftrk,orbit,fmap,ek_ch,re_ch, te_ch) + call tmcat(fsec,re_ch,te_ch,re,te,re,te) + endif + if (fcentre) return call tmdrf(fsec,ftrk,orbit,fmap,dl,ek0,rw,tw) diff --git a/src/util.f90 b/src/util.f90 index 22ed66cfa..8b5f3a193 100644 --- a/src/util.f90 +++ b/src/util.f90 @@ -111,6 +111,7 @@ module code_constfi integer, parameter :: code_nllens = 42 integer, parameter :: code_rfmultipole = 43 integer, parameter :: code_collimator = 44 + integer, parameter :: code_changerefpc = 45 end module code_constfi module aperture_enums From efee29bb952a7bb8c1e8b56ca590222c27d3aac3 Mon Sep 17 00:00:00 2001 From: Tobias Date: Wed, 13 May 2020 13:52:24 +0200 Subject: [PATCH 10/20] small fix. --- src/twiss.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/twiss.f90 b/src/twiss.f90 index 25b8889d0..cbfe60dce 100644 --- a/src/twiss.f90 +++ b/src/twiss.f90 @@ -6513,10 +6513,10 @@ SUBROUTINE tmrf(fsec,ftrk,fcentre,orbit,fmap,el,ds,ek,re,te) call tmcat(fsec,re,te,rw,tw,re,te) - if (changeref_p0) then + call tmchenergy(ftrk,orbit,fmap,ek_ch,re_ch, te_ch) call tmcat(fsec,re_ch,te_ch,re,te,re,te) - endif + if (fcentre) return From a9c0c8418198ac94b6b1e170f238e61c5568acd7 Mon Sep 17 00:00:00 2001 From: Tobias Date: Tue, 26 May 2020 13:54:55 +0200 Subject: [PATCH 11/20] Fix some errors and added for linac. --- libs/DISTlib/source/distgeneration.c | 8 +++- libs/DISTlib/source/distinterface.c | 5 ++- src/mad_dict.c | 13 +++++-- src/mad_dist.c | 56 ++++++++++++++++++++++------ 4 files changed, 64 insertions(+), 18 deletions(-) diff --git a/libs/DISTlib/source/distgeneration.c b/libs/DISTlib/source/distgeneration.c index 544e85c4d..f0fa50da6 100644 --- a/libs/DISTlib/source/distgeneration.c +++ b/libs/DISTlib/source/distgeneration.c @@ -199,7 +199,9 @@ int particle_within_limits_normalized(double *normalized){ } void createcoordinates(int index, double start, double stop, int length, int type){ - double temp [length]; + double *temp; + temp = (double*)malloc(length*sizeof(double)); + if(type ==0){ //Constant value for(int i=0;i incoord[i]->coord[index] = start; @@ -255,6 +257,8 @@ void createcoordinates(int index, double start, double stop, int length, int ty } else issue_error("Unknown type of spacing"); + + free(temp); } @@ -309,6 +313,7 @@ double randray(double mu, double sigma){ void allocateincoord(int linecount){ + printf("allllocating \n"); dist->incoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); dist->outcoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); dist->totincoord = linecount; @@ -330,6 +335,7 @@ void allocateincoord(int linecount){ } dist->isallocated =1; + printf("allllocating2 \n"); } void deallocateincoord(void){ diff --git a/libs/DISTlib/source/distinterface.c b/libs/DISTlib/source/distinterface.c index 95a25f1ff..04ce60530 100644 --- a/libs/DISTlib/source/distinterface.c +++ b/libs/DISTlib/source/distinterface.c @@ -126,7 +126,7 @@ void settasmatrix_element(double value, int row, int column){ void settwisstas(double betax, double alfax, double betay, double alfay){ dist->tas[0][0] = sqrt(betax); - dist->tas[1][0] =-(alfax)/sqrt(betax); + dist->tas[1][0] =-alfax/sqrt(betax); dist->tas[1][1] =-1/sqrt(betax); @@ -167,6 +167,7 @@ void setcoords(double *xn, double *xpn, double *yn, double *ypn, double *zn, dou void settotalsteps(int totgenerate){ dist->totincoord = totgenerate; + printf("setttttinglllll \n"); } void setscan_para_diagonal(int variable, int variable_type, int type, double start, double stop){ @@ -178,7 +179,7 @@ void setscan_para_diagonal(int variable, int variable_type, int type, double sta } dist->ref->typeused[variable] = variable_type; dist->incoordtype=variable_type; // this is wrong - + createcoordinates(variable, start, stop, dist->totincoord,type); } diff --git a/src/mad_dict.c b/src/mad_dict.c index de9963728..c1951ae68 100644 --- a/src/mad_dict.c +++ b/src/mad_dict.c @@ -235,10 +235,15 @@ const char *const_command_def = "horizontal = [s, zero, zero], " "vertical = [s, zero, zero], " "longitudinal = [s, zero, zero], " -"table = [s, disttable, disttable], " -"cutsig_h = [r, {-1}], " -"cutsig_v = [r, {-1}], " -"cutsig_l = [r, {-1}]; " +"table = [s, disttable, disttable], " +"cutsig_h = [r, {-1}], " +"cutsig_v = [r, {-1}], " +"betx = [r, 0], " +"bety = [r, 0], " +"alfx = [r, 0], " +"alfy = [r, 0], " +"use_intial = [l, false, true], " +"cutsig_l = [r, {-1}]; " " " "deselect: control none 0 0 " /* GR: stub exists but not properly implemented */ "flag = [s, none, none], " diff --git a/src/mad_dist.c b/src/mad_dist.c index a706642dc..b1d15055e 100644 --- a/src/mad_dist.c +++ b/src/mad_dist.c @@ -65,19 +65,44 @@ void pro_distribution(struct in_cmd* p){ int cut_l; double cuts[3]; int const C_HOR = 0,C_VER = 2,C_LON = 4; + double *xd, *pxd, *yd, *pyd, *td, *ptd; double npart_d = command_par_value("npart", p->clone); int npart = (int)npart_d; - double x[npart], px[npart], y[npart], py[npart], t[npart], pt[npart]; + xd = mymalloc("distribution", npart * sizeof *xd); + pxd = mymalloc("distribution", npart * sizeof *xd); + yd = mymalloc("distribution", npart * sizeof *xd); + pyd = mymalloc("distribution", npart * sizeof *xd); + td = mymalloc("distribution", npart * sizeof *xd); + ptd = mymalloc("distribution", npart * sizeof *xd); + + + + //double px[npart], y[npart], py[npart], t[npart], pt[npart]; initializedistribution(1); setemitt12(get_value("beam", "ex"), get_value("beam", "ey")); setemitt3(get_value("beam", "et")); sete0andmass0(get_value("beam", "energy"),get_value("beam", "mass") ); - trbegn_(oneturnmat,eigen); - settasmatrixtranspose(eigen); + printf("innnlllll2 \n"); settotalsteps(npart); + printf("innnlllll3 \n"); + if(command_par_value("use_intial", p->clone)!=0){ + double betx = command_par_value("betx", p->clone); + double alfx = command_par_value("alfx", p->clone); + double bety = command_par_value("bety", p->clone); + double alfy = command_par_value("alfy", p->clone); + settwisstas(betx, alfx, bety, alfy); + + printf("ssssss %f %f %f %f \n", betx, alfx, bety, alfy); + } + else{ + trbegn_(oneturnmat,eigen); + settasmatrixtranspose(eigen); + + } + printf("innnlllll4 \n"); type = command_par_string("horizontal", p->clone); cut_l = command_par_vector("cutsig_h", p->clone, cuts); setdistparameters(type, cut_l, cuts, C_HOR, dist_type, start_v,stop_v, coord_t); @@ -89,12 +114,14 @@ void pro_distribution(struct in_cmd* p){ type = command_par_string("longitudinal", p->clone); cut_l = command_par_vector("cutsig_l", p->clone, cuts); setdistparameters(type, cut_l, cuts, C_LON, dist_type, start_v,stop_v, coord_t); - + printf("innnlllll5 \n"); for(int i=0; i<6; i++){ setscan_para_diagonal(i,coord_t[i],dist_type[i],start_v[i],stop_v[i]); } - getunconvertedcoord(x,px,y,py,t,pt, &npart); + + printf("innnlllll6 \n"); + getunconvertedcoord(xd,pxd,yd,pyd,td,ptd, &npart); table_name = command_par_string("table", p->clone); @@ -103,12 +130,12 @@ void pro_distribution(struct in_cmd* p){ for(int i=0; i< npart; i++){ - double_to_table_curr(table_name, "x", &x[i]); - double_to_table_curr(table_name, "px", &px[i]); - double_to_table_curr(table_name, "y", &y[i]); - double_to_table_curr(table_name, "py", &py[i]); - double_to_table_curr(table_name, "t", &t[i]); - double_to_table_curr(table_name, "pt", &pt[i]); + double_to_table_curr(table_name, "x", &xd[i]); + double_to_table_curr(table_name, "px", &pxd[i]); + double_to_table_curr(table_name, "y", &yd[i]); + double_to_table_curr(table_name, "py", &pyd[i]); + double_to_table_curr(table_name, "t", &td[i]); + double_to_table_curr(table_name, "pt", &ptd[i]); num = (int)i; double_to_table_curr(table_name, "number", &num); augmentcountonly(table_name); @@ -136,5 +163,12 @@ void pro_distribution(struct in_cmd* p){ type = command_par_string("longitudinal", p->clone); table_add_header(dist_t, "@ DIST_TYPE_L %%%02ds \"%s\"", strlen(type), stoupper(type)); + //free(xd); + //free(pxd); + //free(yd); + //free(pyd); + //free(td); + //free(ptd); free_distribution(0); // free the memory in the distribution block + } \ No newline at end of file From 39d2fe4012c961852082bddbe34d0616942f24a8 Mon Sep 17 00:00:00 2001 From: Tobias Date: Wed, 27 May 2020 13:45:46 +0200 Subject: [PATCH 12/20] Added changeref0. --- src/trrun.f90 | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/src/trrun.f90 b/src/trrun.f90 index 7526c2eac..a67205593 100644 --- a/src/trrun.f90 +++ b/src/trrun.f90 @@ -1844,8 +1844,37 @@ subroutine ttrf(track,ktrack) enddo endif !! frs add-on end + !call ttchangep0(track,ktrack) end subroutine ttrf +subroutine ttchangep0(track,ktrack) + use math_constfi, only : zero, two, one + use phys_constfi, only : clight + implicit none + double precision :: track(6,*) + double precision :: get_value, bet0 + double precision :: pc0, px_, py_, pt_, onedp + integer :: i, ktrack + + pc0 = get_value('beam ','pc ') + bet0 = get_value('beam ','beta ') + + do i =1, ktrack + px_ = track(1,i) + py_ = track(3,i) + pt_ = track(6,i) + + onedp = sqrt( one + two*pt_/bet0 + (pt_**2)) + + TRACK(2,i) = TRACK(2,i)/onedp + TRACK(4,i) = TRACK(4,i)/onedp + TRACK(6,i) = zero + end do + +end subroutine ttchangep0 + + + subroutine ttcrabrf(track,ktrack,turn) use math_constfi, only : zero, ten3m, ten6p, twopi use phys_constfi, only : clight @@ -1919,7 +1948,8 @@ subroutine ttcrabrf(track,ktrack,turn) TRACK(6,1:ktrack) = TRACK(6,1:ktrack) - & omega * vrf * TRACK(1,1:ktrack) * cos(phirf - bvk*omega*TRACK(5,1:ktrack)) - + + end subroutine ttcrabrf subroutine tthacdip(track,ktrack,turn) From 048fc46361334b3c5679b5a6af74e5c354d461aa Mon Sep 17 00:00:00 2001 From: Tobias Date: Thu, 28 May 2020 15:03:30 +0200 Subject: [PATCH 13/20] Still in some desperate need of cleaning. --- src/madx_ptc_twiss.f90 | 162 +++++++++++++++++++++-------------------- src/trrun.f90 | 2 +- src/twiss.f90 | 1 + 3 files changed, 85 insertions(+), 80 deletions(-) diff --git a/src/madx_ptc_twiss.f90 b/src/madx_ptc_twiss.f90 index b9a28f040..1fd3b81b6 100644 --- a/src/madx_ptc_twiss.f90 +++ b/src/madx_ptc_twiss.f90 @@ -1728,7 +1728,7 @@ subroutine puttwisstable(transfermap,transfermapSaved) double precision :: deltae ! for reference energy increase via acceleration double precision :: deltap ! for deltap treatment ! added on 3 November 2010 to hold Edwards & Teng parametrization - real(dp) :: betx,bety,alfx,alfy,R11,R12,R21,R22 + real(dp) :: betx,bety,alfx,alfy,R11,R12,R21,R22, pt_, onedp ! to convert between Ripken and Edwards-Teng parametrization real(dp) :: kappa,u,ax,ay,kx,ky,kxy2,usqrt,bx,by,cx,cy,cosvp,sinvp,cosvm,sinvm,cosv2,sinv2,cosv1,sinv1 real(dp) :: deltaeValue @@ -1844,76 +1844,80 @@ subroutine puttwisstable(transfermap,transfermapSaved) endif - deltap = A_script_probe%x(5).sub.'0' - deltae = deltae * (1.0 + deltap) - - - opt_fun(beta11)= tw%beta(1,1) * deltae ! beta11=1 - opt_fun(beta12)= tw%beta(1,2) * deltae - opt_fun(beta13)= tw%beta(1,3) * deltae - opt_fun(beta21)= tw%beta(2,1) * deltae - opt_fun(beta22)= tw%beta(2,2) * deltae - opt_fun(beta23)= tw%beta(2,3) * deltae - opt_fun(beta31)= tw%beta(3,1) * deltae - opt_fun(beta32)= tw%beta(3,2) * deltae - opt_fun(beta33)= tw%beta(3,3) * deltae - - opt_fun(alfa11)= tw%alfa(1,1) * deltae - opt_fun(alfa12)= tw%alfa(1,2) * deltae - opt_fun(alfa13)= tw%alfa(1,3) * deltae - opt_fun(alfa21)= tw%alfa(2,1) * deltae - opt_fun(alfa22)= tw%alfa(2,2) * deltae - opt_fun(alfa23)= tw%alfa(2,3) * deltae - opt_fun(alfa31)= tw%alfa(3,1) * deltae - opt_fun(alfa32)= tw%alfa(3,2) * deltae - opt_fun(alfa33)= tw%alfa(3,3) * deltae - - opt_fun(gama11)= tw%gama(1,1) * deltae - opt_fun(gama12)= tw%gama(1,2) * deltae - opt_fun(gama13)= tw%gama(1,3) * deltae - opt_fun(gama21)= tw%gama(2,1) * deltae - opt_fun(gama22)= tw%gama(2,2) * deltae - opt_fun(gama23)= tw%gama(2,3) * deltae - opt_fun(gama31)= tw%gama(3,1) * deltae - opt_fun(gama32)= tw%gama(3,2) * deltae - opt_fun(gama33)= tw%gama(3,3) * deltae + !deltap = A_script_probe%x(5).sub.'0' + !deltae = deltae * (1.0 + deltap) + print *, "deltaeeeee1", deltae, deltap, beta0 + pt_ = A_script_probe%x(5).sub.'0' + onedp = sqrt( one + two*pt_/relativisticBeta + (pt_**2)) + print *, "deltaeeeee", deltae, deltap, relativisticBeta, onedp + + + opt_fun(beta11)= tw%beta(1,1) * onedp ! beta11=1 + opt_fun(beta12)= tw%beta(1,2) * onedp + opt_fun(beta13)= tw%beta(1,3) * onedp + opt_fun(beta21)= tw%beta(2,1) * onedp + opt_fun(beta22)= tw%beta(2,2) * onedp + opt_fun(beta23)= tw%beta(2,3) * onedp + opt_fun(beta31)= tw%beta(3,1) * onedp + opt_fun(beta32)= tw%beta(3,2) * onedp + opt_fun(beta33)= tw%beta(3,3) * onedp + + opt_fun(alfa11)= tw%alfa(1,1) * onedp + opt_fun(alfa12)= tw%alfa(1,2) * onedp + opt_fun(alfa13)= tw%alfa(1,3) * onedp + opt_fun(alfa21)= tw%alfa(2,1) * onedp + opt_fun(alfa22)= tw%alfa(2,2) * onedp + opt_fun(alfa23)= tw%alfa(2,3) * onedp + opt_fun(alfa31)= tw%alfa(3,1) * onedp + opt_fun(alfa32)= tw%alfa(3,2) * onedp + opt_fun(alfa33)= tw%alfa(3,3) * onedp + + opt_fun(gama11)= tw%gama(1,1) * onedp + opt_fun(gama12)= tw%gama(1,2) * onedp + opt_fun(gama13)= tw%gama(1,3) * onedp + opt_fun(gama21)= tw%gama(2,1) * onedp + opt_fun(gama22)= tw%gama(2,2) * onedp + opt_fun(gama23)= tw%gama(2,3) * onedp + opt_fun(gama31)= tw%gama(3,1) * onedp + opt_fun(gama32)= tw%gama(3,2) * onedp + opt_fun(gama33)= tw%gama(3,3) * onedp ! --- derivatives of Twiss paramters w.r.t delta_p - ! NOW why do we need to multiply by deltae, as for the other Twiss parameters? + ! NOW why do we need to multiply by onedp, as for the other Twiss parameters? if (deltap_dependency) then - opt_fun(beta11p)= tw%beta_p(1,1) * deltae - opt_fun(beta12p)= tw%beta_p(1,2) * deltae - opt_fun(beta13p)= tw%beta_p(1,3) * deltae - opt_fun(beta21p)= tw%beta_p(2,1) * deltae - opt_fun(beta22p)= tw%beta_p(2,2) * deltae - opt_fun(beta23p)= tw%beta_p(2,3) * deltae - opt_fun(beta32p)= tw%beta_p(3,2) * deltae - opt_fun(beta33p)= tw%beta_p(3,3) * deltae - - opt_fun(alfa11p)= tw%alfa_p(1,1) * deltae - opt_fun(alfa12p)= tw%alfa_p(1,2) * deltae - opt_fun(alfa13p)= tw%alfa_p(1,3) * deltae - opt_fun(alfa21p)= tw%alfa_p(2,1) * deltae - opt_fun(alfa22p)= tw%alfa_p(2,2) * deltae - opt_fun(alfa23p)= tw%alfa_p(2,3) * deltae - opt_fun(alfa31p)= tw%alfa_p(3,1) * deltae - opt_fun(alfa32p)= tw%alfa_p(3,2) * deltae - opt_fun(alfa33p)= tw%alfa_p(3,3) * deltae - - opt_fun(gama11p)= tw%gama_p(1,1) * deltae - opt_fun(gama12p)= tw%gama_p(1,2) * deltae - opt_fun(gama13p)= tw%gama_p(1,3) * deltae - opt_fun(gama21p)= tw%gama_p(2,1) * deltae - opt_fun(gama22p)= tw%gama_p(2,2) * deltae - opt_fun(gama23p)= tw%gama_p(2,3) * deltae - opt_fun(gama31p)= tw%gama_p(3,1) * deltae - opt_fun(gama32p)= tw%gama_p(3,2) * deltae - opt_fun(gama33p)= tw%gama_p(3,3) * deltae + opt_fun(beta11p)= tw%beta_p(1,1) * onedp + opt_fun(beta12p)= tw%beta_p(1,2) * onedp + opt_fun(beta13p)= tw%beta_p(1,3) * onedp + opt_fun(beta21p)= tw%beta_p(2,1) * onedp + opt_fun(beta22p)= tw%beta_p(2,2) * onedp + opt_fun(beta23p)= tw%beta_p(2,3) * onedp + opt_fun(beta32p)= tw%beta_p(3,2) * onedp + opt_fun(beta33p)= tw%beta_p(3,3) * onedp + + opt_fun(alfa11p)= tw%alfa_p(1,1) * onedp + opt_fun(alfa12p)= tw%alfa_p(1,2) * onedp + opt_fun(alfa13p)= tw%alfa_p(1,3) * onedp + opt_fun(alfa21p)= tw%alfa_p(2,1) * onedp + opt_fun(alfa22p)= tw%alfa_p(2,2) * onedp + opt_fun(alfa23p)= tw%alfa_p(2,3) * onedp + opt_fun(alfa31p)= tw%alfa_p(3,1) * onedp + opt_fun(alfa32p)= tw%alfa_p(3,2) * onedp + opt_fun(alfa33p)= tw%alfa_p(3,3) * onedp + + opt_fun(gama11p)= tw%gama_p(1,1) * onedp + opt_fun(gama12p)= tw%gama_p(1,2) * onedp + opt_fun(gama13p)= tw%gama_p(1,3) * onedp + opt_fun(gama21p)= tw%gama_p(2,1) * onedp + opt_fun(gama22p)= tw%gama_p(2,2) * onedp + opt_fun(gama23p)= tw%gama_p(2,3) * onedp + opt_fun(gama31p)= tw%gama_p(3,1) * onedp + opt_fun(gama32p)= tw%gama_p(3,2) * onedp + opt_fun(gama33p)= tw%gama_p(3,3) * onedp endif ! --- end - ! march 10th: do we need to multiply by deltae the following? + ! march 10th: do we need to multiply by onedp the following? opt_fun(mu1)=tw%mu(1) !* deltae opt_fun(mu2)=tw%mu(2) !* deltae opt_fun(mu3)=tw%mu(3) !* deltae @@ -2034,10 +2038,10 @@ subroutine puttwisstable(transfermap,transfermapSaved) ! to get the same values between twiss and ptc_twiss ! beta11, alfa11 etc... are multiplied by deltae before output ! hence we reflect this in the formula from Lebedev - betx = tw%beta(1,1) * deltaeValue - bety = tw%beta(2,2) * deltaeValue - alfx = tw%alfa(1,1) * deltaeValue - alfy = tw%alfa(2,2) * deltaeValue + betx = tw%beta(1,1) * onedp + bety = tw%beta(2,2) * onedp + alfx = tw%alfa(1,1) * onedp + alfy = tw%alfa(2,2) * onedp else @@ -2045,9 +2049,9 @@ subroutine puttwisstable(transfermap,transfermapSaved) ky=sqrt(tw%beta(2,1)/tw%beta(2,2)); ! beta11, alfa11 etc... are multiplied by deltae before output - ax=kx*tw%alfa(1,1) * deltaeValue -tw%alfa(1,2) * deltaeValue /kx; + ax=kx*tw%alfa(1,1) * onedp -tw%alfa(1,2) * onedp /kx; ! hence we reflect this in the formula from Lebedev - ay=ky*tw%alfa(2,2) * deltaeValue -tw%alfa(2,1) * deltaeValue /ky; + ay=ky*tw%alfa(2,2) * onedp -tw%alfa(2,1) * onedp /ky; kxy2=kx*kx*ky*ky; @@ -2071,10 +2075,10 @@ subroutine puttwisstable(transfermap,transfermapSaved) kappa=one-u - betx = (tw%beta(1,1)/kappa) * deltaeValue - bety = (tw%beta(2,2)/kappa) * deltaeValue - alfx = (tw%alfa(1,1)/kappa) * deltaeValue - alfy = (tw%alfa(2,2)/kappa) * deltaeValue + betx = (tw%beta(1,1)/kappa) * onedp + bety = (tw%beta(2,2)/kappa) * onedp + alfx = (tw%alfa(1,1)/kappa) * onedp + alfy = (tw%alfa(2,2)/kappa) * onedp bx = kx*kappa+u/kx by = ky*kappa-u/ky @@ -2106,10 +2110,10 @@ subroutine puttwisstable(transfermap,transfermapSaved) print*, "ky=",ky print*, "kxy2=",kxy2 - betx = tw%beta(1,1) * deltaeValue - bety = tw%beta(2,2) * deltaeValue - alfx = tw%alfa(1,1) * deltaeValue - alfy = tw%alfa(2,2) * deltaeValue + betx = tw%beta(1,1) * onedp + bety = tw%beta(2,2) * onedp + alfx = tw%alfa(1,1) * onedp + alfy = tw%alfa(2,2) * onedp endif endif @@ -2128,7 +2132,7 @@ subroutine puttwisstable(transfermap,transfermapSaved) ! !--- track the Twiss functions' extremas - call trackBetaExtrema(tw%beta,deltae,betx,bety,tw%disp) + call trackBetaExtrema(tw%beta,onedp,betx,bety,tw%disp) !--- end subroutine puttwisstable diff --git a/src/trrun.f90 b/src/trrun.f90 index 169cc54da..afe952f7f 100644 --- a/src/trrun.f90 +++ b/src/trrun.f90 @@ -1805,7 +1805,7 @@ subroutine ttrf(track,ktrack) enddo endif !! frs add-on end - !call ttchangep0(track,ktrack) + call ttchangep0(track,ktrack) end subroutine ttrf subroutine ttchangep0(track,ktrack) diff --git a/src/twiss.f90 b/src/twiss.f90 index cbfe60dce..6e0e733d6 100644 --- a/src/twiss.f90 +++ b/src/twiss.f90 @@ -6387,6 +6387,7 @@ SUBROUTINE tmchenergy(ftrk,orbit,fmap,ek,re, te) gamma1 = energy1/mass beta1i = energy1/pc1 + print *, "kkkk", pc1/pc re(2,2) = pc/pc1 re(4,4) = pc/pc1 re(6,6) = pc/pc1 From f4da22b316ca152c721491a097468417a8cb14c1 Mon Sep 17 00:00:00 2001 From: tpersson Date: Fri, 29 May 2020 15:53:58 +0200 Subject: [PATCH 14/20] Changes the scaling of the TWISS function from PTC. --- CMakeLists.txt | 0 Makefile_c | 0 libs/DISTlib/CMakeLists.txt | 0 libs/DISTlib/source/CMakeLists.txt | 0 libs/DISTlib/source/distgeneration.c | 0 libs/DISTlib/source/distgeneration.h | 0 libs/DISTlib/source/distinterface.c | 0 libs/DISTlib/source/distinterface.h | 0 libs/DISTlib/source/feigen.c | 0 libs/DISTlib/source/helper.c | 0 libs/DISTlib/source/helper.h | 0 libs/DISTlib/source/outputdist.c | 0 libs/DISTlib/source/outputdist.h | 0 src/CMakeLists.txt | 0 src/mad_cmd.c | 0 src/mad_dict.c | 0 src/mad_dist.c | 0 src/mad_dist.h | 0 src/mad_extrn_f.h | 0 src/mad_gcst.c | 0 src/mad_gcst.h | 0 src/madx.h | 0 src/madx_ptc_twiss.f90 | 91 ++++++++++++++-------------- src/trrun.f90 | 0 src/twiss.f90 | 0 src/util.f90 | 0 tests/test-plot-2/test-plot-2.ref | 0 27 files changed, 45 insertions(+), 46 deletions(-) mode change 100644 => 100755 CMakeLists.txt mode change 100644 => 100755 Makefile_c mode change 100644 => 100755 libs/DISTlib/CMakeLists.txt mode change 100644 => 100755 libs/DISTlib/source/CMakeLists.txt mode change 100644 => 100755 libs/DISTlib/source/distgeneration.c mode change 100644 => 100755 libs/DISTlib/source/distgeneration.h mode change 100644 => 100755 libs/DISTlib/source/distinterface.c mode change 100644 => 100755 libs/DISTlib/source/distinterface.h mode change 100644 => 100755 libs/DISTlib/source/feigen.c mode change 100644 => 100755 libs/DISTlib/source/helper.c mode change 100644 => 100755 libs/DISTlib/source/helper.h mode change 100644 => 100755 libs/DISTlib/source/outputdist.c mode change 100644 => 100755 libs/DISTlib/source/outputdist.h mode change 100644 => 100755 src/CMakeLists.txt mode change 100644 => 100755 src/mad_cmd.c mode change 100644 => 100755 src/mad_dict.c mode change 100644 => 100755 src/mad_dist.c mode change 100644 => 100755 src/mad_dist.h mode change 100644 => 100755 src/mad_extrn_f.h mode change 100644 => 100755 src/mad_gcst.c mode change 100644 => 100755 src/mad_gcst.h mode change 100644 => 100755 src/madx.h mode change 100644 => 100755 src/madx_ptc_twiss.f90 mode change 100644 => 100755 src/trrun.f90 mode change 100644 => 100755 src/twiss.f90 mode change 100644 => 100755 src/util.f90 mode change 100644 => 100755 tests/test-plot-2/test-plot-2.ref diff --git a/CMakeLists.txt b/CMakeLists.txt old mode 100644 new mode 100755 diff --git a/Makefile_c b/Makefile_c old mode 100644 new mode 100755 diff --git a/libs/DISTlib/CMakeLists.txt b/libs/DISTlib/CMakeLists.txt old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/CMakeLists.txt b/libs/DISTlib/source/CMakeLists.txt old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/distgeneration.c b/libs/DISTlib/source/distgeneration.c old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/distgeneration.h b/libs/DISTlib/source/distgeneration.h old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/distinterface.c b/libs/DISTlib/source/distinterface.c old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/distinterface.h b/libs/DISTlib/source/distinterface.h old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/feigen.c b/libs/DISTlib/source/feigen.c old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/helper.c b/libs/DISTlib/source/helper.c old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/helper.h b/libs/DISTlib/source/helper.h old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/outputdist.c b/libs/DISTlib/source/outputdist.c old mode 100644 new mode 100755 diff --git a/libs/DISTlib/source/outputdist.h b/libs/DISTlib/source/outputdist.h old mode 100644 new mode 100755 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt old mode 100644 new mode 100755 diff --git a/src/mad_cmd.c b/src/mad_cmd.c old mode 100644 new mode 100755 diff --git a/src/mad_dict.c b/src/mad_dict.c old mode 100644 new mode 100755 diff --git a/src/mad_dist.c b/src/mad_dist.c old mode 100644 new mode 100755 diff --git a/src/mad_dist.h b/src/mad_dist.h old mode 100644 new mode 100755 diff --git a/src/mad_extrn_f.h b/src/mad_extrn_f.h old mode 100644 new mode 100755 diff --git a/src/mad_gcst.c b/src/mad_gcst.c old mode 100644 new mode 100755 diff --git a/src/mad_gcst.h b/src/mad_gcst.h old mode 100644 new mode 100755 diff --git a/src/madx.h b/src/madx.h old mode 100644 new mode 100755 diff --git a/src/madx_ptc_twiss.f90 b/src/madx_ptc_twiss.f90 old mode 100644 new mode 100755 index 1fd3b81b6..9275b221b --- a/src/madx_ptc_twiss.f90 +++ b/src/madx_ptc_twiss.f90 @@ -1846,10 +1846,9 @@ subroutine puttwisstable(transfermap,transfermapSaved) !deltap = A_script_probe%x(5).sub.'0' !deltae = deltae * (1.0 + deltap) - print *, "deltaeeeee1", deltae, deltap, beta0 pt_ = A_script_probe%x(5).sub.'0' onedp = sqrt( one + two*pt_/relativisticBeta + (pt_**2)) - print *, "deltaeeeee", deltae, deltap, relativisticBeta, onedp + opt_fun(beta11)= tw%beta(1,1) * onedp ! beta11=1 @@ -1862,25 +1861,25 @@ subroutine puttwisstable(transfermap,transfermapSaved) opt_fun(beta32)= tw%beta(3,2) * onedp opt_fun(beta33)= tw%beta(3,3) * onedp - opt_fun(alfa11)= tw%alfa(1,1) * onedp - opt_fun(alfa12)= tw%alfa(1,2) * onedp - opt_fun(alfa13)= tw%alfa(1,3) * onedp - opt_fun(alfa21)= tw%alfa(2,1) * onedp - opt_fun(alfa22)= tw%alfa(2,2) * onedp - opt_fun(alfa23)= tw%alfa(2,3) * onedp - opt_fun(alfa31)= tw%alfa(3,1) * onedp - opt_fun(alfa32)= tw%alfa(3,2) * onedp - opt_fun(alfa33)= tw%alfa(3,3) * onedp - - opt_fun(gama11)= tw%gama(1,1) * onedp - opt_fun(gama12)= tw%gama(1,2) * onedp - opt_fun(gama13)= tw%gama(1,3) * onedp - opt_fun(gama21)= tw%gama(2,1) * onedp - opt_fun(gama22)= tw%gama(2,2) * onedp - opt_fun(gama23)= tw%gama(2,3) * onedp - opt_fun(gama31)= tw%gama(3,1) * onedp - opt_fun(gama32)= tw%gama(3,2) * onedp - opt_fun(gama33)= tw%gama(3,3) * onedp + opt_fun(alfa11)= tw%alfa(1,1) + opt_fun(alfa12)= tw%alfa(1,2) + opt_fun(alfa13)= tw%alfa(1,3) + opt_fun(alfa21)= tw%alfa(2,1) + opt_fun(alfa22)= tw%alfa(2,2) + opt_fun(alfa23)= tw%alfa(2,3) + opt_fun(alfa31)= tw%alfa(3,1) + opt_fun(alfa32)= tw%alfa(3,2) + opt_fun(alfa33)= tw%alfa(3,3) + + opt_fun(gama11)= tw%gama(1,1) / onedp + opt_fun(gama12)= tw%gama(1,2) / onedp + opt_fun(gama13)= tw%gama(1,3) / onedp + opt_fun(gama21)= tw%gama(2,1) / onedp + opt_fun(gama22)= tw%gama(2,2) / onedp + opt_fun(gama23)= tw%gama(2,3) / onedp + opt_fun(gama31)= tw%gama(3,1) / onedp + opt_fun(gama32)= tw%gama(3,2) / onedp + opt_fun(gama33)= tw%gama(3,3) / onedp ! --- derivatives of Twiss paramters w.r.t delta_p @@ -1895,25 +1894,25 @@ subroutine puttwisstable(transfermap,transfermapSaved) opt_fun(beta32p)= tw%beta_p(3,2) * onedp opt_fun(beta33p)= tw%beta_p(3,3) * onedp - opt_fun(alfa11p)= tw%alfa_p(1,1) * onedp - opt_fun(alfa12p)= tw%alfa_p(1,2) * onedp - opt_fun(alfa13p)= tw%alfa_p(1,3) * onedp - opt_fun(alfa21p)= tw%alfa_p(2,1) * onedp - opt_fun(alfa22p)= tw%alfa_p(2,2) * onedp - opt_fun(alfa23p)= tw%alfa_p(2,3) * onedp - opt_fun(alfa31p)= tw%alfa_p(3,1) * onedp - opt_fun(alfa32p)= tw%alfa_p(3,2) * onedp - opt_fun(alfa33p)= tw%alfa_p(3,3) * onedp - - opt_fun(gama11p)= tw%gama_p(1,1) * onedp - opt_fun(gama12p)= tw%gama_p(1,2) * onedp - opt_fun(gama13p)= tw%gama_p(1,3) * onedp - opt_fun(gama21p)= tw%gama_p(2,1) * onedp - opt_fun(gama22p)= tw%gama_p(2,2) * onedp - opt_fun(gama23p)= tw%gama_p(2,3) * onedp - opt_fun(gama31p)= tw%gama_p(3,1) * onedp - opt_fun(gama32p)= tw%gama_p(3,2) * onedp - opt_fun(gama33p)= tw%gama_p(3,3) * onedp + opt_fun(alfa11p)= tw%alfa_p(1,1) + opt_fun(alfa12p)= tw%alfa_p(1,2) + opt_fun(alfa13p)= tw%alfa_p(1,3) + opt_fun(alfa21p)= tw%alfa_p(2,1) + opt_fun(alfa22p)= tw%alfa_p(2,2) + opt_fun(alfa23p)= tw%alfa_p(2,3) + opt_fun(alfa31p)= tw%alfa_p(3,1) + opt_fun(alfa32p)= tw%alfa_p(3,2) + opt_fun(alfa33p)= tw%alfa_p(3,3) + + opt_fun(gama11p)= tw%gama_p(1,1) / onedp + opt_fun(gama12p)= tw%gama_p(1,2) / onedp + opt_fun(gama13p)= tw%gama_p(1,3) / onedp + opt_fun(gama21p)= tw%gama_p(2,1) / onedp + opt_fun(gama22p)= tw%gama_p(2,2) / onedp + opt_fun(gama23p)= tw%gama_p(2,3) / onedp + opt_fun(gama31p)= tw%gama_p(3,1) / onedp + opt_fun(gama32p)= tw%gama_p(3,2) / onedp + opt_fun(gama33p)= tw%gama_p(3,3) / onedp endif ! --- end @@ -2040,8 +2039,8 @@ subroutine puttwisstable(transfermap,transfermapSaved) ! hence we reflect this in the formula from Lebedev betx = tw%beta(1,1) * onedp bety = tw%beta(2,2) * onedp - alfx = tw%alfa(1,1) * onedp - alfy = tw%alfa(2,2) * onedp + alfx = tw%alfa(1,1) + alfy = tw%alfa(2,2) else @@ -2077,8 +2076,8 @@ subroutine puttwisstable(transfermap,transfermapSaved) betx = (tw%beta(1,1)/kappa) * onedp bety = (tw%beta(2,2)/kappa) * onedp - alfx = (tw%alfa(1,1)/kappa) * onedp - alfy = (tw%alfa(2,2)/kappa) * onedp + alfx = (tw%alfa(1,1)/kappa) + alfy = (tw%alfa(2,2)/kappa) bx = kx*kappa+u/kx by = ky*kappa-u/ky @@ -2112,8 +2111,8 @@ subroutine puttwisstable(transfermap,transfermapSaved) betx = tw%beta(1,1) * onedp bety = tw%beta(2,2) * onedp - alfx = tw%alfa(1,1) * onedp - alfy = tw%alfa(2,2) * onedp + alfx = tw%alfa(1,1) + alfy = tw%alfa(2,2) endif endif diff --git a/src/trrun.f90 b/src/trrun.f90 old mode 100644 new mode 100755 diff --git a/src/twiss.f90 b/src/twiss.f90 old mode 100644 new mode 100755 diff --git a/src/util.f90 b/src/util.f90 old mode 100644 new mode 100755 diff --git a/tests/test-plot-2/test-plot-2.ref b/tests/test-plot-2/test-plot-2.ref old mode 100644 new mode 100755 From af842de6ee2e99afcc74b86b430963bf31c8ff58 Mon Sep 17 00:00:00 2001 From: tpersson Date: Fri, 5 Jun 2020 11:20:20 +0200 Subject: [PATCH 15/20] Some cleanup. --- src/mad_dist.c | 9 +++------ src/mad_option.c | 4 ++-- src/trrun.f90 | 4 ++-- src/twiss.f90 | 7 ++----- 4 files changed, 9 insertions(+), 15 deletions(-) diff --git a/src/mad_dist.c b/src/mad_dist.c index b1d15055e..257132a1b 100755 --- a/src/mad_dist.c +++ b/src/mad_dist.c @@ -84,9 +84,7 @@ void pro_distribution(struct in_cmd* p){ setemitt12(get_value("beam", "ex"), get_value("beam", "ey")); setemitt3(get_value("beam", "et")); sete0andmass0(get_value("beam", "energy"),get_value("beam", "mass") ); - printf("innnlllll2 \n"); settotalsteps(npart); - printf("innnlllll3 \n"); if(command_par_value("use_intial", p->clone)!=0){ double betx = command_par_value("betx", p->clone); @@ -95,14 +93,13 @@ void pro_distribution(struct in_cmd* p){ double alfy = command_par_value("alfy", p->clone); settwisstas(betx, alfx, bety, alfy); - printf("ssssss %f %f %f %f \n", betx, alfx, bety, alfy); } else{ trbegn_(oneturnmat,eigen); settasmatrixtranspose(eigen); } - printf("innnlllll4 \n"); + type = command_par_string("horizontal", p->clone); cut_l = command_par_vector("cutsig_h", p->clone, cuts); setdistparameters(type, cut_l, cuts, C_HOR, dist_type, start_v,stop_v, coord_t); @@ -114,13 +111,13 @@ void pro_distribution(struct in_cmd* p){ type = command_par_string("longitudinal", p->clone); cut_l = command_par_vector("cutsig_l", p->clone, cuts); setdistparameters(type, cut_l, cuts, C_LON, dist_type, start_v,stop_v, coord_t); - printf("innnlllll5 \n"); + for(int i=0; i<6; i++){ setscan_para_diagonal(i,coord_t[i],dist_type[i],start_v[i],stop_v[i]); } - printf("innnlllll6 \n"); + getunconvertedcoord(xd,pxd,yd,pyd,td,ptd, &npart); diff --git a/src/mad_option.c b/src/mad_option.c index b589e6ca8..ebc103eb3 100644 --- a/src/mad_option.c +++ b/src/mad_option.c @@ -1,6 +1,6 @@ #include "madx.h" -int + int get_option(const char* str) { /* This function is called by fortran to get option of a command */ @@ -40,7 +40,7 @@ set_defaults(const char* str) /* reset options, beam etc. to defaults */ if ((pos = name_list_pos(str, defined_commands->list)) > -1) { if (strcmp(str, "option") == 0) - { + { if (options != NULL) delete_command(options); options = clone_command(defined_commands->commands[pos]); } diff --git a/src/trrun.f90 b/src/trrun.f90 index afe952f7f..f8ab0e1ed 100755 --- a/src/trrun.f90 +++ b/src/trrun.f90 @@ -1061,7 +1061,7 @@ SUBROUTINE ttmult_cf_mini(track,ktrack,dxt,dyt,turn, thin_foc) dipr = bvk * normal(0) !vals(1,0) dipi = bvk * skew(0) - !!print *, "ggggg_old", gstr, sstr + ! cf magnet with quadrupole & sextupole gstr = normal(1)/elrad sstr = normal(2)/elrad @@ -1805,7 +1805,7 @@ subroutine ttrf(track,ktrack) enddo endif !! frs add-on end - call ttchangep0(track,ktrack) + end subroutine ttrf subroutine ttchangep0(track,ktrack) diff --git a/src/twiss.f90 b/src/twiss.f90 index 6e0e733d6..87731aee9 100755 --- a/src/twiss.f90 +++ b/src/twiss.f90 @@ -3443,7 +3443,6 @@ SUBROUTINE tmmap(code,fsec,ftrk,orbit,fmap,ek,re,te,fcentre,dl) ! nothing on purpose! case (code_wire) - call tmchenergy(ftrk,orbit,fmap,ek,re, te) ! nothing for now... case (code_dipedge) @@ -4627,10 +4626,10 @@ SUBROUTINE tmmult(fsec,ftrk,orbit,fmap,re,te) di = (dr * y + di * x) / (iord) + f_errors(2*iord+1) dr = drt enddo - ! re(2,1) = - dr + re(2,1) = - dr re(2,3) = + di re(4,1) = + di - ! re(4,3) = + dr + re(4,3) = + dr endif !---- Add the missing focussing component of thin dipoles @@ -4673,7 +4672,6 @@ SUBROUTINE tmmult(fsec,ftrk,orbit,fmap,re,te) te(4,3,3) = - di endif endif -te = zero end SUBROUTINE tmmult SUBROUTINE tmoct(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te) @@ -6387,7 +6385,6 @@ SUBROUTINE tmchenergy(ftrk,orbit,fmap,ek,re, te) gamma1 = energy1/mass beta1i = energy1/pc1 - print *, "kkkk", pc1/pc re(2,2) = pc/pc1 re(4,4) = pc/pc1 re(6,6) = pc/pc1 From ef43891e4e9c33a0bf6d995da34560d0e9bb2fdf Mon Sep 17 00:00:00 2001 From: tpersson Date: Fri, 5 Jun 2020 15:22:03 +0200 Subject: [PATCH 16/20] Small changes. --- .travis.yml | 5 ++++- Makefile_test | 1 + tests/test-match-10/line_out.txt.cfg | 4 ++-- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/.travis.yml b/.travis.yml index aa933fb60..995b58de0 100644 --- a/.travis.yml +++ b/.travis.yml @@ -15,7 +15,10 @@ install: - make numdiff-linux64-gnu - pip install cython - - git clone https://github.com/hibtc/cpymad.git + - git clone https://github.com/tpersson/cpymad.git + - cd cpymad + - git checkout addingDistribution + - cd .. - ./buildcpymad.sh - cd cpymad && python setup.py build_ext --no-X11 --blas --lapack --madxdir=$MADXDIR - pip install -e . diff --git a/Makefile_test b/Makefile_test index ac2928bd2..a9cc0c5b7 100644 --- a/Makefile_test +++ b/Makefile_test @@ -47,6 +47,7 @@ test-memory test-control \ test-makethin test-makethin-2 test-makethin-3 test-makethin-4 test-elseparator\ test-survey test-survey-2 test-survey-3 \ test-track test-track-2 test-track-3 test-track-4 test-track-5 test-track-6 \ +test-track-7 test-track-8 test-track-9 test-track-10 test-track-11 test-track-12 test-track-13 test-track-14 \ test-track-acd test-track-rotations \ test-twiss test-twiss-2 test-twiss-3 test-twiss-4 test-twiss-5 \ test-twiss-6 test-twiss-8 test-twiss-9 test-twiss-10 test-twiss-11 \ diff --git a/tests/test-match-10/line_out.txt.cfg b/tests/test-match-10/line_out.txt.cfg index 38ea1055a..bb25ded58 100644 --- a/tests/test-match-10/line_out.txt.cfg +++ b/tests/test-match-10/line_out.txt.cfg @@ -1,4 +1,4 @@ 44-46 * skip # head -24-26 1 abs=2e-9 # dq1, dq2 -30-32 1 abs=1e-8 # betx, bety +24-26 1 abs=4e-9 # dq1, dq2 +30-32 1 abs=2e-8 # betx, bety 48-$ * any abs=6e-10 rel=4e-7 From 16dfaf692b38eee614fdade1c242c437b947a7e8 Mon Sep 17 00:00:00 2001 From: tpersson Date: Fri, 5 Jun 2020 15:35:31 +0200 Subject: [PATCH 17/20] Make these changes in a separate request. --- src/madx_ptc_twiss.f90 | 161 ++++++++++++++++++++--------------------- 1 file changed, 79 insertions(+), 82 deletions(-) diff --git a/src/madx_ptc_twiss.f90 b/src/madx_ptc_twiss.f90 index 9275b221b..b9a28f040 100755 --- a/src/madx_ptc_twiss.f90 +++ b/src/madx_ptc_twiss.f90 @@ -1728,7 +1728,7 @@ subroutine puttwisstable(transfermap,transfermapSaved) double precision :: deltae ! for reference energy increase via acceleration double precision :: deltap ! for deltap treatment ! added on 3 November 2010 to hold Edwards & Teng parametrization - real(dp) :: betx,bety,alfx,alfy,R11,R12,R21,R22, pt_, onedp + real(dp) :: betx,bety,alfx,alfy,R11,R12,R21,R22 ! to convert between Ripken and Edwards-Teng parametrization real(dp) :: kappa,u,ax,ay,kx,ky,kxy2,usqrt,bx,by,cx,cy,cosvp,sinvp,cosvm,sinvm,cosv2,sinv2,cosv1,sinv1 real(dp) :: deltaeValue @@ -1844,79 +1844,76 @@ subroutine puttwisstable(transfermap,transfermapSaved) endif - !deltap = A_script_probe%x(5).sub.'0' - !deltae = deltae * (1.0 + deltap) - pt_ = A_script_probe%x(5).sub.'0' - onedp = sqrt( one + two*pt_/relativisticBeta + (pt_**2)) - - - - opt_fun(beta11)= tw%beta(1,1) * onedp ! beta11=1 - opt_fun(beta12)= tw%beta(1,2) * onedp - opt_fun(beta13)= tw%beta(1,3) * onedp - opt_fun(beta21)= tw%beta(2,1) * onedp - opt_fun(beta22)= tw%beta(2,2) * onedp - opt_fun(beta23)= tw%beta(2,3) * onedp - opt_fun(beta31)= tw%beta(3,1) * onedp - opt_fun(beta32)= tw%beta(3,2) * onedp - opt_fun(beta33)= tw%beta(3,3) * onedp - - opt_fun(alfa11)= tw%alfa(1,1) - opt_fun(alfa12)= tw%alfa(1,2) - opt_fun(alfa13)= tw%alfa(1,3) - opt_fun(alfa21)= tw%alfa(2,1) - opt_fun(alfa22)= tw%alfa(2,2) - opt_fun(alfa23)= tw%alfa(2,3) - opt_fun(alfa31)= tw%alfa(3,1) - opt_fun(alfa32)= tw%alfa(3,2) - opt_fun(alfa33)= tw%alfa(3,3) - - opt_fun(gama11)= tw%gama(1,1) / onedp - opt_fun(gama12)= tw%gama(1,2) / onedp - opt_fun(gama13)= tw%gama(1,3) / onedp - opt_fun(gama21)= tw%gama(2,1) / onedp - opt_fun(gama22)= tw%gama(2,2) / onedp - opt_fun(gama23)= tw%gama(2,3) / onedp - opt_fun(gama31)= tw%gama(3,1) / onedp - opt_fun(gama32)= tw%gama(3,2) / onedp - opt_fun(gama33)= tw%gama(3,3) / onedp + deltap = A_script_probe%x(5).sub.'0' + deltae = deltae * (1.0 + deltap) + + + opt_fun(beta11)= tw%beta(1,1) * deltae ! beta11=1 + opt_fun(beta12)= tw%beta(1,2) * deltae + opt_fun(beta13)= tw%beta(1,3) * deltae + opt_fun(beta21)= tw%beta(2,1) * deltae + opt_fun(beta22)= tw%beta(2,2) * deltae + opt_fun(beta23)= tw%beta(2,3) * deltae + opt_fun(beta31)= tw%beta(3,1) * deltae + opt_fun(beta32)= tw%beta(3,2) * deltae + opt_fun(beta33)= tw%beta(3,3) * deltae + + opt_fun(alfa11)= tw%alfa(1,1) * deltae + opt_fun(alfa12)= tw%alfa(1,2) * deltae + opt_fun(alfa13)= tw%alfa(1,3) * deltae + opt_fun(alfa21)= tw%alfa(2,1) * deltae + opt_fun(alfa22)= tw%alfa(2,2) * deltae + opt_fun(alfa23)= tw%alfa(2,3) * deltae + opt_fun(alfa31)= tw%alfa(3,1) * deltae + opt_fun(alfa32)= tw%alfa(3,2) * deltae + opt_fun(alfa33)= tw%alfa(3,3) * deltae + + opt_fun(gama11)= tw%gama(1,1) * deltae + opt_fun(gama12)= tw%gama(1,2) * deltae + opt_fun(gama13)= tw%gama(1,3) * deltae + opt_fun(gama21)= tw%gama(2,1) * deltae + opt_fun(gama22)= tw%gama(2,2) * deltae + opt_fun(gama23)= tw%gama(2,3) * deltae + opt_fun(gama31)= tw%gama(3,1) * deltae + opt_fun(gama32)= tw%gama(3,2) * deltae + opt_fun(gama33)= tw%gama(3,3) * deltae ! --- derivatives of Twiss paramters w.r.t delta_p - ! NOW why do we need to multiply by onedp, as for the other Twiss parameters? + ! NOW why do we need to multiply by deltae, as for the other Twiss parameters? if (deltap_dependency) then - opt_fun(beta11p)= tw%beta_p(1,1) * onedp - opt_fun(beta12p)= tw%beta_p(1,2) * onedp - opt_fun(beta13p)= tw%beta_p(1,3) * onedp - opt_fun(beta21p)= tw%beta_p(2,1) * onedp - opt_fun(beta22p)= tw%beta_p(2,2) * onedp - opt_fun(beta23p)= tw%beta_p(2,3) * onedp - opt_fun(beta32p)= tw%beta_p(3,2) * onedp - opt_fun(beta33p)= tw%beta_p(3,3) * onedp - - opt_fun(alfa11p)= tw%alfa_p(1,1) - opt_fun(alfa12p)= tw%alfa_p(1,2) - opt_fun(alfa13p)= tw%alfa_p(1,3) - opt_fun(alfa21p)= tw%alfa_p(2,1) - opt_fun(alfa22p)= tw%alfa_p(2,2) - opt_fun(alfa23p)= tw%alfa_p(2,3) - opt_fun(alfa31p)= tw%alfa_p(3,1) - opt_fun(alfa32p)= tw%alfa_p(3,2) - opt_fun(alfa33p)= tw%alfa_p(3,3) - - opt_fun(gama11p)= tw%gama_p(1,1) / onedp - opt_fun(gama12p)= tw%gama_p(1,2) / onedp - opt_fun(gama13p)= tw%gama_p(1,3) / onedp - opt_fun(gama21p)= tw%gama_p(2,1) / onedp - opt_fun(gama22p)= tw%gama_p(2,2) / onedp - opt_fun(gama23p)= tw%gama_p(2,3) / onedp - opt_fun(gama31p)= tw%gama_p(3,1) / onedp - opt_fun(gama32p)= tw%gama_p(3,2) / onedp - opt_fun(gama33p)= tw%gama_p(3,3) / onedp + opt_fun(beta11p)= tw%beta_p(1,1) * deltae + opt_fun(beta12p)= tw%beta_p(1,2) * deltae + opt_fun(beta13p)= tw%beta_p(1,3) * deltae + opt_fun(beta21p)= tw%beta_p(2,1) * deltae + opt_fun(beta22p)= tw%beta_p(2,2) * deltae + opt_fun(beta23p)= tw%beta_p(2,3) * deltae + opt_fun(beta32p)= tw%beta_p(3,2) * deltae + opt_fun(beta33p)= tw%beta_p(3,3) * deltae + + opt_fun(alfa11p)= tw%alfa_p(1,1) * deltae + opt_fun(alfa12p)= tw%alfa_p(1,2) * deltae + opt_fun(alfa13p)= tw%alfa_p(1,3) * deltae + opt_fun(alfa21p)= tw%alfa_p(2,1) * deltae + opt_fun(alfa22p)= tw%alfa_p(2,2) * deltae + opt_fun(alfa23p)= tw%alfa_p(2,3) * deltae + opt_fun(alfa31p)= tw%alfa_p(3,1) * deltae + opt_fun(alfa32p)= tw%alfa_p(3,2) * deltae + opt_fun(alfa33p)= tw%alfa_p(3,3) * deltae + + opt_fun(gama11p)= tw%gama_p(1,1) * deltae + opt_fun(gama12p)= tw%gama_p(1,2) * deltae + opt_fun(gama13p)= tw%gama_p(1,3) * deltae + opt_fun(gama21p)= tw%gama_p(2,1) * deltae + opt_fun(gama22p)= tw%gama_p(2,2) * deltae + opt_fun(gama23p)= tw%gama_p(2,3) * deltae + opt_fun(gama31p)= tw%gama_p(3,1) * deltae + opt_fun(gama32p)= tw%gama_p(3,2) * deltae + opt_fun(gama33p)= tw%gama_p(3,3) * deltae endif ! --- end - ! march 10th: do we need to multiply by onedp the following? + ! march 10th: do we need to multiply by deltae the following? opt_fun(mu1)=tw%mu(1) !* deltae opt_fun(mu2)=tw%mu(2) !* deltae opt_fun(mu3)=tw%mu(3) !* deltae @@ -2037,10 +2034,10 @@ subroutine puttwisstable(transfermap,transfermapSaved) ! to get the same values between twiss and ptc_twiss ! beta11, alfa11 etc... are multiplied by deltae before output ! hence we reflect this in the formula from Lebedev - betx = tw%beta(1,1) * onedp - bety = tw%beta(2,2) * onedp - alfx = tw%alfa(1,1) - alfy = tw%alfa(2,2) + betx = tw%beta(1,1) * deltaeValue + bety = tw%beta(2,2) * deltaeValue + alfx = tw%alfa(1,1) * deltaeValue + alfy = tw%alfa(2,2) * deltaeValue else @@ -2048,9 +2045,9 @@ subroutine puttwisstable(transfermap,transfermapSaved) ky=sqrt(tw%beta(2,1)/tw%beta(2,2)); ! beta11, alfa11 etc... are multiplied by deltae before output - ax=kx*tw%alfa(1,1) * onedp -tw%alfa(1,2) * onedp /kx; + ax=kx*tw%alfa(1,1) * deltaeValue -tw%alfa(1,2) * deltaeValue /kx; ! hence we reflect this in the formula from Lebedev - ay=ky*tw%alfa(2,2) * onedp -tw%alfa(2,1) * onedp /ky; + ay=ky*tw%alfa(2,2) * deltaeValue -tw%alfa(2,1) * deltaeValue /ky; kxy2=kx*kx*ky*ky; @@ -2074,10 +2071,10 @@ subroutine puttwisstable(transfermap,transfermapSaved) kappa=one-u - betx = (tw%beta(1,1)/kappa) * onedp - bety = (tw%beta(2,2)/kappa) * onedp - alfx = (tw%alfa(1,1)/kappa) - alfy = (tw%alfa(2,2)/kappa) + betx = (tw%beta(1,1)/kappa) * deltaeValue + bety = (tw%beta(2,2)/kappa) * deltaeValue + alfx = (tw%alfa(1,1)/kappa) * deltaeValue + alfy = (tw%alfa(2,2)/kappa) * deltaeValue bx = kx*kappa+u/kx by = ky*kappa-u/ky @@ -2109,10 +2106,10 @@ subroutine puttwisstable(transfermap,transfermapSaved) print*, "ky=",ky print*, "kxy2=",kxy2 - betx = tw%beta(1,1) * onedp - bety = tw%beta(2,2) * onedp - alfx = tw%alfa(1,1) - alfy = tw%alfa(2,2) + betx = tw%beta(1,1) * deltaeValue + bety = tw%beta(2,2) * deltaeValue + alfx = tw%alfa(1,1) * deltaeValue + alfy = tw%alfa(2,2) * deltaeValue endif endif @@ -2131,7 +2128,7 @@ subroutine puttwisstable(transfermap,transfermapSaved) ! !--- track the Twiss functions' extremas - call trackBetaExtrema(tw%beta,onedp,betx,bety,tw%disp) + call trackBetaExtrema(tw%beta,deltae,betx,bety,tw%disp) !--- end subroutine puttwisstable From e98ab8e197b1f1b41f958b4decc8ab9572d7c8f3 Mon Sep 17 00:00:00 2001 From: tpersson Date: Fri, 5 Jun 2020 15:57:48 +0200 Subject: [PATCH 18/20] Removed a test in RF-cavity --- src/twiss.f90 | 4 ++-- tests/test-thick-quad-2/test-thick-quad-2.madx | 6 +++++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/twiss.f90 b/src/twiss.f90 index cf7234e67..e1787964c 100755 --- a/src/twiss.f90 +++ b/src/twiss.f90 @@ -6490,8 +6490,8 @@ SUBROUTINE tmrf(fsec,ftrk,fcentre,orbit,fmap,el,ds,ek,re,te) call tmcat(fsec,re,te,rw,tw,re,te) - call tmchenergy(ftrk,orbit,fmap,ek_ch,re_ch, te_ch) - call tmcat(fsec,re_ch,te_ch,re,te,re,te) + ! call tmchenergy(ftrk,orbit,fmap,ek_ch,re_ch, te_ch) + !call tmcat(fsec,re_ch,te_ch,re,te,re,te) if (fcentre) return diff --git a/tests/test-thick-quad-2/test-thick-quad-2.madx b/tests/test-thick-quad-2/test-thick-quad-2.madx index 306b03ab8..ffeb2e72f 100644 --- a/tests/test-thick-quad-2/test-thick-quad-2.madx +++ b/tests/test-thick-quad-2/test-thick-quad-2.madx @@ -39,13 +39,17 @@ rflag:= 1/2; rflag:=3/8; tr$macro(turn): macro = { +value, rfvolt; } track, onepass, update, dump; start, x=1e-6, px=0, y=1e-6, py=0, t=-0.25; start, x=1e-6, px=0, y=1e-6, py=0, t=-1.0; start, x=1e-6, px=0, y=1e-6, py=0, t=-1.75; -run, turns=50, ffile=1; +start, x=1e-6, px=0, y=1e-6, py=0, t=-0.25; +start, x=1e-6, px=0, y=1e-6, py=0, t=-1.0; +start, x=1e-6, px=0, y=1e-6, py=0, t=-1.75; +run, turns=150, ffile=1; endtrack; stop; From 05e7050b098ac2a871660cdd5ddbf8af12ce8a2e Mon Sep 17 00:00:00 2001 From: tpersson Date: Fri, 5 Jun 2020 16:08:31 +0200 Subject: [PATCH 19/20] Restore the original thick quad test. --- tests/test-thick-quad-2/test-thick-quad-2.madx | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/tests/test-thick-quad-2/test-thick-quad-2.madx b/tests/test-thick-quad-2/test-thick-quad-2.madx index ffeb2e72f..306b03ab8 100644 --- a/tests/test-thick-quad-2/test-thick-quad-2.madx +++ b/tests/test-thick-quad-2/test-thick-quad-2.madx @@ -39,17 +39,13 @@ rflag:= 1/2; rflag:=3/8; tr$macro(turn): macro = { -value, rfvolt; } track, onepass, update, dump; start, x=1e-6, px=0, y=1e-6, py=0, t=-0.25; start, x=1e-6, px=0, y=1e-6, py=0, t=-1.0; start, x=1e-6, px=0, y=1e-6, py=0, t=-1.75; -start, x=1e-6, px=0, y=1e-6, py=0, t=-0.25; -start, x=1e-6, px=0, y=1e-6, py=0, t=-1.0; -start, x=1e-6, px=0, y=1e-6, py=0, t=-1.75; -run, turns=150, ffile=1; +run, turns=50, ffile=1; endtrack; stop; From 7f106f1f3a6b1f340567f7c73bac11f929f1c027 Mon Sep 17 00:00:00 2001 From: tpersson Date: Mon, 8 Jun 2020 08:53:26 +0200 Subject: [PATCH 20/20] Small cleanup. --- libs/DISTlib/source/distgeneration.c | 5 +---- src/twiss.f90 | 10 +++++----- src/util.f90 | 2 +- 3 files changed, 7 insertions(+), 10 deletions(-) diff --git a/libs/DISTlib/source/distgeneration.c b/libs/DISTlib/source/distgeneration.c index f0fa50da6..6bb3cc336 100755 --- a/libs/DISTlib/source/distgeneration.c +++ b/libs/DISTlib/source/distgeneration.c @@ -252,7 +252,6 @@ void createcoordinates(int index, double start, double stop, int length, int ty tmp = randray(start, stop); dist->incoord[i]->coord[index] = tmp; - //printf("%f \n", sqrt(tmp/2)); } } else @@ -313,7 +312,6 @@ double randray(double mu, double sigma){ void allocateincoord(int linecount){ - printf("allllocating \n"); dist->incoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); dist->outcoord = (struct coordinates**)malloc(linecount*sizeof(struct coordinates*)); dist->totincoord = linecount; @@ -335,7 +333,6 @@ void allocateincoord(int linecount){ } dist->isallocated =1; - printf("allllocating2 \n"); } void deallocateincoord(void){ @@ -343,6 +340,6 @@ void deallocateincoord(void){ free(dist->incoord); free(dist->outcoord); //dist->ref->typeused = (int*)malloc(dim*sizeof(int)); - dist->isallocated =0; + dist->isallocated = 0; dist->totincoord = -1; } \ No newline at end of file diff --git a/src/twiss.f90 b/src/twiss.f90 index e1787964c..56fd21802 100755 --- a/src/twiss.f90 +++ b/src/twiss.f90 @@ -3465,8 +3465,8 @@ SUBROUTINE tmmap(code,fsec,ftrk,orbit,fmap,ek,re,te,fcentre,dl) case (code_rfmultipole) call tmrfmult(fsec,ftrk,orbit,fmap,ek,re,te) - case (code_changerefpc) - call tmchenergy(ftrk,orbit,fmap,ek,re, te) + case (code_changerefp0) + call tmchp0(ftrk,orbit,fmap,ek,re, te) case default !--- anything else: ! nil (23, 28, 34) @@ -6321,7 +6321,7 @@ SUBROUTINE tmdrf(fsec,ftrk,orbit,fmap,dl,ek,re,te) end SUBROUTINE tmdrf -SUBROUTINE tmchenergy(ftrk,orbit,fmap,ek,re, te) +SUBROUTINE tmchp0(ftrk,orbit,fmap,ek,re, te) use twisslfi use twiss_elpfi use twissbeamfi, only : deltap, pc, gamma, energy, beta @@ -6372,7 +6372,7 @@ SUBROUTINE tmchenergy(ftrk,orbit,fmap,ek,re, te) if(ftrk) call tmtrak(ek,re,te,orbit,orbit) -end SUBROUTINE tmchenergy +end SUBROUTINE tmchp0 SUBROUTINE tmrf(fsec,ftrk,fcentre,orbit,fmap,el,ds,ek,re,te) @@ -6490,7 +6490,7 @@ SUBROUTINE tmrf(fsec,ftrk,fcentre,orbit,fmap,el,ds,ek,re,te) call tmcat(fsec,re,te,rw,tw,re,te) - ! call tmchenergy(ftrk,orbit,fmap,ek_ch,re_ch, te_ch) + ! call tmchp0(ftrk,orbit,fmap,ek_ch,re_ch, te_ch) !call tmcat(fsec,re_ch,te_ch,re,te,re,te) diff --git a/src/util.f90 b/src/util.f90 index 8b5f3a193..141de9370 100755 --- a/src/util.f90 +++ b/src/util.f90 @@ -111,7 +111,7 @@ module code_constfi integer, parameter :: code_nllens = 42 integer, parameter :: code_rfmultipole = 43 integer, parameter :: code_collimator = 44 - integer, parameter :: code_changerefpc = 45 + integer, parameter :: code_changerefp0 = 45 end module code_constfi module aperture_enums