Changeset 14306


Ignore:
Timestamp:
08/02/12 13:22:38 (12 years ago)
Author:
tbretz
Message:
Fixed a performance issue with the calculation of the moon position, convert the light condition directly to currents, renamed the light-condition to current-prediction.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/FACT++/src/smartfact.cc

    r14305 r14306  
    18551855
    18561856#ifdef HAVE_NOVA
     1857
     1858    vector<pair<ln_equ_posn, double>> fMoonCoords;
     1859
     1860    void CalcMoonCoords(double jd)
     1861    {
     1862        jd = floor(jd);
     1863
     1864        fMoonCoords.clear();
     1865        for (double h=0; h<1; h+=1./(24*12))
     1866        {
     1867            ln_equ_posn  moon;
     1868            ln_get_lunar_equ_coords_prec(jd+h, &moon, 0.01);
     1869
     1870            const double disk = ln_get_lunar_disk(jd+h);
     1871
     1872            fMoonCoords.push_back(make_pair(moon, disk));
     1873        }
     1874    }
     1875
    18571876    pair<vector<float>, pair<Time, float>> GetVisibility(ln_equ_posn *src, ln_lnlat_posn *observer, double jd)
    18581877    {
     
    18701889        int cnt = 0;
    18711890
     1891        int i=0;
     1892
    18721893        vector<float> alt;
    1873         for (double h=0; h<1; h+=1./(24*12))
     1894        for (double h=0; h<1; h+=1./(24*12), i++)
    18741895        {
    18751896            if (src==0)
    1876                 ln_get_lunar_equ_coords(jd+h, &moon);
     1897                moon = fMoonCoords[i].first; //ln_get_lunar_equ_coords_prec(jd+h, &moon, 0.01);
    18771898
    18781899            ln_hrz_posn hrz;
     
    19051926        const double jd1 = fmod(fSun.fRiseAstronomical.JD(), 1);
    19061927
    1907         ln_equ_posn  moon;
    1908 
    19091928        double max   = -1;
    19101929        double maxjd =  0;
    19111930
    19121931        int cnt = 0;
    1913 
    1914         vector<float> alt;
    1915         for (double h=0; h<1; h+=1./(24*12))
    1916         {
    1917             double lc = -1;
     1932        int i = 0;
     1933
     1934        vector<float> vec;
     1935        for (double h=0; h<1; h+=1./(24*12), i++)
     1936        {
     1937            double cur = -1;
    19181938
    19191939            if (h>jd0 && h<jd1)
    19201940            {
    1921                 const double disk = ln_get_lunar_disk(jd+h);
    1922                 ln_get_lunar_equ_coords(jd+h, &moon);
     1941                ln_equ_posn moon  = fMoonCoords[i].first;//ln_get_lunar_equ_coords_prec(jd+h, &moon, 0.01);
     1942                const double disk = fMoonCoords[i].second;//ln_get_lunar_disk(jd+h);
    19231943
    19241944                ln_hrz_posn hrz;
     
    19311951                const double angle = m.Angle(src->ra, src->dec);
    19321952
    1933                 lc = angle*hrz.alt*disk*M_PI*M_PI/180/180;
    1934 
    1935                 alt.push_back(lc);
     1953                const double lc = angle*hrz.alt*pow(disk, 6)/360/360;
     1954
     1955                cur = 7.7+4942*lc;
     1956
     1957                vec.push_back(cur); // Covert LC to pixel current in uA
    19361958            }
    19371959
    1938             if (lc>max)
     1960            if (cur>max)
    19391961            {
    1940                 max   = lc;
     1962                max   = cur;
    19411963                maxjd = jd+h;
    19421964            }
    19431965
    1944             if (h>jd0 && h<jd1 && lc>0)
     1966            if (h>jd0 && h<jd1 && cur>0)
    19451967                cnt++;
    19461968        }
    19471969
    1948         if (max<=15 || cnt==0)
     1970        if (max<0 || cnt==0)
    19491971            return make_pair(vector<float>(), make_pair(Time(), 0));
    19501972
    1951         return make_pair(alt, make_pair(maxjd, maxjd>jd+jd0&&maxjd<jd+jd1?max:0));
     1973        return make_pair(vec, make_pair(maxjd, maxjd>jd+jd0&&maxjd<jd+jd1?max:0));
    19521974    }
    19531975#endif
     
    19591981
    19601982        Time now;
     1983
     1984        CalcMoonCoords(now.JD());
    19611985
    19621986        fSun  = Sun (lon, lat, now);
     
    20202044        map<Time, pair<string, float>> lightcond;
    20212045        vector<vector<float>> alt;
    2022         vector<vector<float>> hlc;
     2046        vector<vector<float>> cur;
    20232047
    20242048#ifdef HAVE_NOVA
     
    20692093                if (lc.first.size()>0)
    20702094                {
    2071                     hlc.push_back(lc.first);
     2095                    cur.push_back(lc.first);
    20722096                    lightcond[lc.second.first] = make_pair(name, lc.second.second);
    20732097                }
     
    21112135            }
    21122136
     2137            out4 << setprecision(3);
    21132138            for (auto it=lightcond.begin(); it!=lightcond.end(); it++)
    21142139            {
     
    21172142                out4 << "<B>" << it->second.first << "</B>";
    21182143                if (it->second.second>0)
    2119                     out4 << " [" << nearbyint(90-it->second.second) << "&deg;]";
     2144                    out4 << " [" << nearbyint(it->second.second) << "]";
    21202145            }
    21212146
     
    21342159            out2 << HTML::kWhite << '\t' << Time()-now << '\n';
    21352160
    2136             WriteBinaryVec(now, "hist-visibility",      alt, 75, 15, "Alt "+title.str());
    2137             WriteBinaryVec(now, "hist-light-condition", hlc,  1,  0, "LC "+title.str());
     2161            WriteBinaryVec(now, "hist-visibility",         alt, 75, 15, "Alt "+title.str());
     2162            WriteBinaryVec(now, "hist-current-prediction", cur, 100,  0, "I "  +title.str());
    21382163        }
    21392164        catch (const exception &e)
     
    21512176        ofstream(fPath+"/source-list.data") << out2.str();
    21522177        ofstream(fPath+"/visibility.data") << out3.str();
    2153         ofstream(fPath+"/light-condition.data") << out4.str();
     2178        ofstream(fPath+"/current-prediction.data") << out4.str();
    21542179    }
    21552180
Note: See TracChangeset for help on using the changeset viewer.