Commit 9d1e9a42eb42d6d018f4a4a19656f440d58cc6e7

  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Committer)
  • Sun Jul 10 18:43:22 EEST 2011
  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Author)
  • Sun Jul 10 18:43:22 EEST 2011
Use log(n) bisection-like method to find ktsqr and yindexes using
one single routine.
src/amplitude.cpp
(37 / 21)
  
9696
9797REAL Amplitude::N(REAL ktsqr, REAL y, bool bspline, bool derivative)
9898{
99 if (ktsqr*0.9999999 > MaxKtsqr() or ktsqr*1.00000001 < MinKtsqr())
99 REAL lnktsqr = std::log(ktsqr);
100 if (lnktsqr > lnktsqrvals[KtsqrPoints()-1] or lnktsqr < lnktsqrvals[0])
100101 cerr << "Ktsqrval " << ktsqr << " is out of range [" <<
101102 MinKtsqr() << ", " << MaxKtsqr() << "]! " << LINEINFO << endl;
102103
111111
112112 // Find ktsqrval and yval indexes refer to index for which
113113 // val[index] is smaller than (or equal) y/ktsqr
114 int yind = -1;
115 int ktsqrind = KtsqrIndex(ktsqr);
114 int ktsqrind = FindIndex(ktsqr, ktsqrvals);
115 int yind = FindIndex(y, yvals);
116
116117 for (unsigned int i=0; i<yvals.size()-1; i++)
117118 {
118119 if (yvals[i]<=y and yvals[i+1]>y)
122122 break;
123123 }
124124 }
125 if (yind < 0) // Didn't find, so refers to the largest one
125 if (yind < 0 or ktsqrind<0)
126126 {
127 yind=yvals.size()-1;
128 if (y-yvals[yind]>0.05)
129 cerr << "Asked amplitude at too large Y=" << y << ", falling back to "
130 << " y=" << yvals[yind] << ". " << LINEINFO << endl;
127 cerr << "Something very weird happened at " << LINEINFO << " with "
128 << "ktsqr=" << ktsqr << ", y=" << y << endl;
129 return 0;
131130 }
132131
133132 // Keep y fixed, interpolate ktsqr
155155 // with that
156156 ///TODO: Why?
157157 if (ktsqrind>0 and interpolation_start==0) { interpolation_start=1; }
158
159 //cout << "interp range " << interpolation_start << " - " << interpolation_end << " index " << ktsqrind << endl;
160158
161159 int interpo_points = interpolation_end - interpolation_start+1;
162160
257257REAL Amplitude::LogLogDerivative(REAL ktsqr, REAL rapidity)
258258{
259259
260 int ktsqrind = KtsqrIndex(ktsqr);
260 int ktsqrind = FindIndex(ktsqr, ktsqrvals);
261261 if (ktsqrind==0 or ktsqrind >= KtsqrPoints())
262262 {
263263 cerr << "Can't calculate derivative at the edge of the ktsqr range. "
644644unsigned int Amplitude::YPoints()
645645{
646646 return yvals.size()-1;
647 //return static_cast<unsigned int>(maxy/delta_y)+1;
648647}
649648
650649unsigned int Amplitude::KtsqrPoints()
718718}
719719
720720/* Returns index i for which
721 * ktsqrvals[i]<=ktsqr<ktsqrvals[i+1]
721 * vec[i]<=val
722722 * If such index can't be found, returns -1
723723 */
724int Amplitude::KtsqrIndex(REAL ktsqr)
724
725int Amplitude::FindIndex(REAL val, std::vector<REAL> &vec)
725726{
726 int ktsqrind=-1;
727 if (ktsqr< MinKtsqr() or ktsqr > MaxKtsqr()) return -1;
728 for (unsigned int i=0; i<KtsqrPoints()-1; i++)
727 if (val < vec[0]) return -1;
728
729 int ind=-1;
730
731 uint start=0; uint end=vec.size()-1;
732 while(end-start>10)
729733 {
730 if (ktsqrvals[i]<=ktsqr and ktsqrvals[i+1]>ktsqr)
734 int tmp = static_cast<int>((start+end)/2.0);
735
736 if (vec[tmp]>=val)
737 end=tmp;
738 else
739 start=tmp;
740 }
741
742
743 for (uint i=start; i<=end; i++)
744 {
745 if (vec[i]<=val and vec[i+1]>val)
731746 {
732 ktsqrind=i;
747 ind=i;
733748 break;
734749 }
735750 }
736 if (ktsqrind==-1) return KtsqrPoints()-1;
737 return ktsqrind;
751 if (ind == -1) return vec.size()-1;
752 return ind;
738753}
754
755
756
739757
740758/*
741759 * Add new rapidity value for tables
src/amplitude.hpp
(1 / 1)
  
9494 virtual REAL MinKtsqr();
9595 virtual REAL MaxKtsqr();
9696
97 int KtsqrIndex(REAL ktsqr);
97 int FindIndex(REAL val, std::vector<REAL> &vec);
9898
9999
100100 REAL InitialCondition(REAL ktsqr); // N() at y=0
src/main.cpp
(2 / 2)
  
560560 cout << "# ktsqr amplitude initial_condition bspline_amplitude" << endl;
561561 ktsqr_mult = std::pow(N->KtsqrMultiplier(),1.0/4.0);
562562
563 for (int i=0; i<N->KtsqrPoints()-1; i+=10)
563 for (int i=0; i<N->KtsqrPoints()-1; i++)
564564 {
565 REAL tmpktsqr = N->MinKtsqr()*std::pow(ktsqr_mult, i);
565 REAL tmpktsqr = N->Ktsqrval(i);
566566 cout << std::scientific << std::setprecision(15) << tmpktsqr << " " <<
567567 N->N(tmpktsqr, y) << " " << N->InitialCondition(tmpktsqr)
568568 << " " << N->N(tmpktsqr, y, true) << endl;
src/solver_force.cpp
(1 / 1)
  
470470 */
471471REAL BruteForceSolver::InterpolateN(REAL ktsqr, const REAL* array, bool bspline, bool der)
472472{
473 int ktsqrind = KtsqrIndex(ktsqr);
473 int ktsqrind = FindIndex(ktsqr, ktsqrvals);
474474 if (ktsqr >= MaxKtsqr()) // Didn't find, so refers to the largest one
475475 return array[ktsqrvals.size()-1];
476476 if (ktsqr <= MinKtsqr())