#include <iostream>
#include <TCanvas.h>
#include <TF2.h>
#include "TModel.h"
#include "TBitmap.h"
#include "showtime.C"

TCanvas* c2;
TF2* img2;

TModel* refa = new TModel(1);
TModel* sig1a = new TModel(1);
TModel* sig2a = new TModel(1);
TModel* sig1ph = new TModel(2);
TModel* sig2ph = new TModel(2);
TBitmap* map=new TBitmap(50,50);
TBitmap* data=new TBitmap(*map);

void crank()
{
   UInt_t width=data->GetWidth();
   UInt_t height=data->GetHeight();
   TASImage im("data.jpg");
   im.Scale(width,height);
   TArrayD* ar=im.GetArray();
   for (UInt_t i=0,k=0; i<height; i++) {
      for (UInt_t j=0; j<width; j++,k++) {
         (*data)[i][j] = (*ar)[k];
      }
   }
   delete ar;
   refa->Randomize(10);
   sig1a->Randomize(4);
   *sig2a = *sig1a;
   *sig2a *= 0.83;
   sig1ph->Randomize(30);
   sig2ph->Randomize(30);
   map->Precompute(5);
   map->Interfere(*refa,*sig1a,*sig2a,*sig1ph,*sig2ph);
   map->Draw();
}

Double_t crankmap2(Double_t* var, Double_t* par)
{
   int k=0;
   for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++,k++) {
         if (par[k] == 1) {
            Double_t* p=(*refa)[i];
            p[j] = var[0];
         }
         else if (par[k] == 2) {
            Double_t* p=(*refa)[i];
            p[j] = var[1];
         }
         else {
            Double_t* p=(*refa)[i];
            p[j] = par[k];
         }
      }
   }
   for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++,k++) {
         if (par[k] == 1) {
            Double_t* p=(*sig1a)[i];
            p[j] = var[0];
         }
         else if (par[k] == 2) {
            Double_t* p=(*sig1a)[i];
            p[j] = var[1];
         }
         else {
            Double_t* p=(*sig1a)[i];
            p[j] = par[k];
         }
      }
   }
   for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++,k++) {
         if (par[k] == 1) {
            Double_t* p=(*sig2a)[i];
            p[j] = var[0];
         }
         else if (par[k] == 2) {
            Double_t* p=(*sig2a)[i];
            p[j] = var[1];
         }
         else {
            Double_t* p=(*sig2a)[i];
            p[j] = par[k];
         }
      }
   }
   for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++,k++) {
         if (par[k] == 1) {
            Double_t* p=(*sig1ph)[i];
            p[j] = var[0];
         }
         else if (par[k] == 2) {
            Double_t* p=(*sig1ph)[i];
            p[j] = var[1];
         }
         else {
            Double_t* p=(*sig1ph)[i];
            p[j] = par[k];
         }
      }
   }
   for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++,k++) {
         if (par[k] == 1) {
            Double_t* p=(*sig2ph)[i];
            p[j] = var[0];
         }
         else if (par[k] == 2) {
            Double_t* p=(*sig2ph)[i];
            p[j] = var[1];
         }
         else {
            Double_t* p=(*sig2ph)[i];
            p[j] = par[k];
         }
      }
   }
   map->Interfere(*refa,*sig1a,*sig2a,*sig1ph,*sig2ph);
   return map->ChiSquared(*data);
}

Int_t drawTMap(Int_t par1, Int_t par2, Double_t lowlim1, Double_t highlim1,
                                       Double_t lowlim2, Double_t highlim2)
{
   if (c2) {
      delete c2;
   }
   c2 = new TCanvas("c2","graph window",200,10,700,500);
   if (img2) {
      delete img2;
   }
   img2 = new TF2("img2",crankmap2,lowlim1,highlim1,lowlim2,highlim2,125);
   int k=0;
   for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++,k++) {
         img2->SetParameter(k,(*refa)[i][j]);
         img2->SetParameter(k+25,(*sig1a)[i][j]);
         img2->SetParameter(k+50,(*sig2a)[i][j]);
         img2->SetParameter(k+75,(*sig1ph)[i][j]);
         img2->SetParameter(k+100,(*sig2ph)[i][j]);
      }
   }
   img2->SetParameter(par1,1);
   img2->SetParameter(par2,2);
   img2->SetNpx(10);
   img2->SetNpy(10);
   img2->Draw("surf1");
   c2->Update();
   return 0;
}