!! These three procedures should be adjusted to the concrete situation !! !!!!!!!!!!!!!!! MAIN PROGRAM S P I D E R !!!!!!!!!!!!!!!!!!!!!!!!!!! !! The program calls Spider Subroutines: SPIDRSUB.FOR or SPIDER.SUB !! !! The program uses files: SPIDAT (data) and SPIDER.IN (parameters) !! Cc 22 Nov 2008 Cc I will modify the code created in June-July 2005 to study trapping Cc of the electrons based on comparison between BATSE HXR data and NoRH Cc spatially resolved light curves. Cc The idea is to use the results of numerical solution of the Fokker-Plank Cc Equation (provided by V.F.Melnikov) as input fo the phenomenological Cc treatment in terms of characteristic life-time \tau(E,mu,S) (Fleishman, 2005, Nobe Proc). Cc ! We do not consider emission here, only fast electron evolution. Cc The GOAL is to derive the cube of the life time as a function of E,mu, and S, Cc find an analytical approximation for it, and then use these results for faster modeling Cc of the fast electron evolution in the magnetic traps. This later be implemented into Cc GS simulator to be developed by Gelu Nita and I. Cc First step is to implement integration over irregular time grid inherited from Cc Victor's computations to describe trapping. Start from a "blind" approach: Cc the injection function has gaussian shape with unknown parameters. Cc 25 Nov 2008 Cc The code works appropriately, when: 1) the injection core is explicitly specified, i.e., Cc two-parameter fitting is used (blind [4- and 3- par] test is unstable, perhapse because Cc the injection profile is not too different from the trapping response); 2) the tail of the Cc time profiles must be cut at the level above second time point at the rise phase; Cc otherwise the time profiles can contain a long tail at low level (I do not know yet its origin). Cc The accuracy of fits varies from fit to fit: this means that I must save the parameter errors Cc for further fitting of them by model functions. Cc 07 Dec 2008 Cc The input has been changed to the array (30,60,60,25) with a constant time resolution of 4 sec. Cc 25 time steps. Now I plan to perform global test of the fit with the output life-time(29,60,60) Cc and Ampl(29,60,60) [E(30) returns zeros] Cc 05 Jan 2009 Cc This is a next step fitting process to replace the numeric solution of the F-P equation Cc by the life-time. This code reads the array of the life-time(Mu,S,E) and fit the dependence on Cc Mu by a model function with a few free parameters functions of S & E Cc 17 June 2005 Cc Start to modify the fitting program to find the charachteristic Cc life time of the radio emission as observed along the spatially Cc resolved (in 2d) radio source. Plan: use X-ray light curves as Cc the proxi of the injection profile and fit radio profile to the Cc injection profile convolved with the trapping function Cc First step: test - one - four light curves produced by calculation Cc within trapping model (Anis. paper, Astro-ph/0505145, 2005): gaussian Cc injection profile. Cc Current version uses calculated radio flux at the foot-point (fourth Cc column of the file 'trap_model.dat') as the injection proxi. Fitting is made Cc for one localized light curve (one of four columns of the above file) Cc Q: why the convolution gives relatively poor fit for J=1-3? Although Cc excellent fit for J=4 (self-fit) Cc 20 June 2005 Cc Upgrade of Version 1. 20 local light curves produced by model GS_trap_model. Cc Expected: life time vs position along the loop. z=0.1 - loop top, Cc z=2 - foot point. We'll study spatial dependence of the life time along the loop Cc Accross: thin loop (one pixel). Cc 06 July 2005 Cc Start test the code with original data: X-ray BATSE and 17 GHz Nobeyama Cc (provided by Tim Bastian). The data points (X and R) are taken at different Cc times and have different cadence: 1 s for radio and 1.024 s for X-rays. Cc This version of the code must take this into accont and determine the life time Cc of the radio emission coming from a few different spatial locations. cc We will modify primarily fitting function and reading. Cc 07 July 2005 Cc Start to modify to calculate 2D distributions of the life time and response intensity Cc We first modify all required arrays for the case of 7 Oct 1992 flare, e.g., Cc AllData(43,40,358) Cc 12 July 2005 Cc Previous version is the final working code for the event 7 Oct 1992 (Stox V fit). Cc Now we change the code for 27 Oct 1992, including AllDAta(43x40x358), whic is enough Cc for all four sample bursts Cc 14 July 2005 Cc This version is adjusted from 'Oct_27' to the event 23 Dec 93 Cc 21 July 2005 Cc This version is adjusted from 'Jun_27' back to the event 07 Oct 92: Cc more numbers obtained directy from the IDL 'read-write' code Cc 01 April 2019 Cc This is a dll version of the code. The most recent change is the ability of the code to Cc accept any combination of polarizations: I, L&R, I&V, I&P Cc NO EXPLICIT provision is made for XX or YY only. If true XX(YY) is provided, it muct be supploied twice (as L and R inputs) Cc If input is I then only the Fitting Mode I is allowed Cc If a combination of two polarization inputs is provided, the user can select a different fitting mode. real*8 function Get_Tables(argc, argv) !(ParmIn,Spectrum,ar3) implicit none !GCC$ ATTRIBUTES DLLEXPORT :: get_tables integer, intent(in) :: argc !integer*8, intent(in) :: argc integer, intent(in), dimension(argc) :: argv !integer*8, intent(in), dimension(argc) :: argv !External forstr !call Call_Parms(%val(argv(1)), %val(argv(2))) !call forstr(%val(argv(1)), %val(argv(2))) Get_Tables=2. Open(8,File='Long_input.txt') !Write(8,126)'parm=', ' Default value', ' Name, Units' Write(8,77) 'Nparms;', 7, '; ;user ;Number of fit Parms' Write(8,77) 'Angular Code;', 0, '; ;data ;0L for PK, 1L for FK' Write(8,77) 'Npix;', 128, '; ;data ;Number of pixels sent to dll' !Write(8,77) 'NpicY;', 128, ';Number of pixels along y axes' Write(8,77)'Nfreq;', 24, '; ;data ;Number of frequencies in the spectrum' Write(8,77)'Fitting Mode;', 1, '; ;user ;Case of the fit: I:1, L&R:2, I&V:3, I&P:4' Write(8,77)'Stokes Data;', 1, '; ;data ;Case of the data: I:1, L&R:2, I&V:3, I&P:4' close(8) Open(8,File='Real_input.txt') !Write(8,126)'parm=', ' Default value', ' Name, Units' Write(8,7) 'SIMPLEX Step;', 0.17, '; ;user', ';SIMPLEX Step' Write(8,7) 'SIMPLEX EPS;', 1d-6, '; ;user', ';SIMPLEX accuracy' Write(8,7)'Flux Threshold;',1.,';sfu*GHz',' ; user ;Flux threshold to be fitted' Write(8,7)'Pixel Area;',4.,';arcsec^2',' ;data ;Number of pixels along y axes' Write(8,7)'LOS Depth;',8.,';arcsec ',';user ;LOS depth (assumed)' Write(8,7)'E_min;',0.015d0,';MeV ;user',';Min energy in PLW (assumed)' close(8) Open(8,File='Parms_input.txt') !Write(8,126)'parm=', ' Default value', ' Name, Units' Write(8,8)'n_nth;',1d0,';',1.d-4,';',2d3,';','1d7 cm^-3' &, ';Non-thermal density' Write(8,8) 'B;', 4d0,';',1.d-2,';',3d1,';', '1d2G ' &, ';Magnetic field' Write(8,8)'theta;',60.,';',22.d0,';',87d0,';','deg', &';Viewing angle to B' Write(8,8)'n_th;',10d0,';',1d-2,';',6d2,';', '1d9 cm^-3', & ';Thermal density' Write(8,8)'Delta;',4.5,';',1.7d0,';',15d0,';','No', & ';PLW spectral index' Write(8,8)'E_max;',5d0,';',0.1d0,';',1d1,';','MeV', &';Max energy in PLW' Write(8,8) 'T_e;',5d0,';',1.5d0,';',6d1,';','MK',';Temperature' Write(8,8)'Res1;',5d0,';',0.2,';',0.2d2,';', ' ',';Not used' Write(8,8)'Res2;',5d0,';',0.2,';',0.2d2,';', ' ',';Not used' Write(8,8)'Res3;',5d0,';',0.2,';',0.2d2,';', ' ',';Not used' Write(8,8)'Res4;',5d0,';',0.2,';',0.2d2,';', ' ',';Not used' Write(8,8)'Res5;',5d0,';',0.2,';',0.2d2,';', ' ',';Not used' Write(8,8)'Res6;',5d0,';',0.2,';',0.2d2,';', ' ',';Not used' Write(8,8)'Res7;',5d0,';',0.2,';',0.2d2,';', ' ',';Not used' Write(8,8)'Res8;',5d0,';',0.2,';',0.2d2,';', ' ',';Not used' close(8) Open(8,File='Parms_out.txt') !Write(8,126)'parm=', ' Default value', ' Name, Units' Write(8,7)'n_nth;',1d7,';cm^-3' &, ';Non-thermal density \pm err' Write(8,7) 'B;', 4d2, ';G ' &, ';Magnetic field \pm err' Write(8,7)'theta;',60.,';deg', & ';Viewing angle to B \pm err' Write(8,7)'n_th;',1d10, ';cm^-3', & ';Thermal density \pm err' Write(8,7)'Delta;',4.5,'; ', & ';PLW spectral index \pm err' Write(8,7)'E_max;',5d0,';MeV', &';Max energy in PLW \pm err' Write(8,7) 'T_e;',5d0,';MK',';Temperature \pm err' Write(8,7)'Chi-2;',1d0, '; ',';Chi-2 or Average Fit Error (Dim.less)' close(8) c126 FORMAT(1x,a10,2x,a20,2x,a10,a30) 7 FORMAT(1x,a15,2x,G10.3,2x,a10,a40) 8 FORMAT(1x,a11,3(2x,G10.3,a1),2x,a10,a40) 77 FORMAT(1x,a13,2x,G10.3,2x,a40) !character*8 string(3,29) !string(1,1)='area' end real*8 function GET_MW_FIT(argc, argv) implicit none interface subroutine SPIDER(N_dat_in,r_inp,ParmsIn,fr_in,spec_in,ArrParms,ErrParms,spec_out) INTEGER*8 :: N_dat_in,r_inp,ParmsIn,fr_in,spec_in,ArrParms,ErrParms,spec_out end end interface INTEGER*8:: argc, argv(*) !! !integer, intent(in) :: argc !! integer*8, intent(in) :: argc ! for x64 !! !integer, intent(in), dimension(argc) :: argv !! integer*8, intent(in), dimension(argc) :: argv ! real*8 status,test8,cpa0 !!integer, intent(in) :: arg2 !open(9,file='D:\GS_Modeling\DLL_Tested\Test.dat') !write(9,*)argc !integer, intent(in), dimension(argc) :: argvv !call klein4dll_sub(%val(argv(1)), %val(argv(2)), %val(argv(3))) !call Get_Parms(%val(argv(1)), %val(argv(2))) C common /status/ status C common /test8/ test8 !returns Nthreads for test purpose C common /testMu/cpa0 !returns mu_0 for test purpose call SPIDER (%val(argv(1)), & %val(argv(2)), %val(argv(3)), %val(argv(4)), %val(argv(5)) & , %val(argv(6)), %val(argv(7)), %val(argv(8))) get_mw_fit=0d0 C transfer_sol_Slice(%val(argv(1)), %val(argv(2)), C & %val(argv(3)), %val(argv(4))) c call transfer_sol(%val(argv(1)), %val(argv(2)), c & %val(argv(3))) !Get_MW_Slice= cpa0 !status !1 return end !program subroutine & SPIDER(N_dat_in,r_inp,ParmsIn,fr_in,spec_in,ArrParms,ErrParms &,spec_out) !(xxx_in,yyy_in) ! Fitting of Func(X1,X2). IMPLICIT REAL*8 (A-H,O-Z) character KEY common/VAR/EPS,STEP,NPAR,NX1,NX2,KX common/ARR/FUNDAT(300,9),ERRDAT(300,9),X1dat(300),X2dat(9), ! NX1max=300, NX2max=9 *PAR(15),DPAR(15),VECT(15),PSCALE(15),PARmin(15),PARmax(15), ! NPARmax=15 *FUNC(19),AS(19,15) ! 19=NPARmax+4 (for SIMPLEX) c common/AllDat/AllData(43,40,358) common/Indecs/Npic,Nfreq common/Eminmax/E_min, E_max, T_fix,are_arc,dep_arc common/B2Ind/NB2 common /NType/IAType common/Fit_case/N_case_d,N_case_f C common/ParGuess/ !Integer*4 N_dat_in(4) Integer*4 N_dat_in(6) Real*8 ParGuess(15),r_inp(6),fr_in(N_dat_in(4)),ParmsIn(15,3), & ArrParms(N_dat_in(3),8), & ErrParms(N_dat_in(3),8), & spec_in(N_dat_in(3),N_dat_in(4),4), & spec_out(N_dat_in(3),N_dat_in(4),2), & FreqData(N_dat_in(4)),Flcp(N_dat_in(3),N_dat_in(4)), & Frcp(N_dat_in(3),N_dat_in(4)), & ErrFlcp(N_dat_in(3),N_dat_in(4)),ErrFrcp(N_dat_in(3),N_dat_in(4)) & ,FFitR(N_dat_in(3),N_dat_in(4)),FFitL(N_dat_in(3),N_dat_in(4)) &,ArrNrl(N_dat_in(3)),ArrBgauss(N_dat_in(3)), & ArrEmax(N_dat_in(3)),ArrAngle(N_dat_in(3)) &,ArrNth(N_dat_in(3)),ArrDelta1(N_dat_in(3)) &,ErrNrl(N_dat_in(3)),ErrBgauss(N_dat_in(3)), & ErrEmax(N_dat_in(3)),ErrAngle(N_dat_in(3)) &,ErrNth(N_dat_in(3)),ErrDelta1(N_dat_in(3)),ErrData(N_dat_in(3)), & RChi2(N_dat_in(3)) &,ArrT(N_dat_in(3)),ErrT(N_dat_in(3)) C & FFitR(N_dat_in(3),N_dat_in(4),N_dat_in(5)), C & FFitL(N_dat_in(3),N_dat_in(4),N_dat_in(5)), C C Real*8, ALLOCATABLE :: FreqData(:),Flcp(:,:),Frcp(:,:), C &ErrFlcp(:,:),ErrFrcp(:,:),FFitR(:,:),FFitL(:,:) C &,ArrNrl(:),ArrBgauss(:),ArrEmax(:),ArrAngle(:) C &,ArrNth(:),ArrDelta1(:) C &,ErrNrl(:),ErrBgauss(:),ErrEmax(:),ErrAngle(:) C &,ErrNth(:),ErrDelta1(:),ErrData(:),RChi2(:) C &,ArrT(:),ErrT(:) CCC &,ArrParms(:,:,:),ErrParms(:,:,:) C Real*4 ParMu1(30,60),ParMu2(30,60),ParMu3(30,60) C &,ErrParMu1(30,60),ErrParMu2(30,60),ErrParMu3(30,60) ! increase dimentions of the arrays later for the 2D case - done on 07 Dec 2008 NX41=60 NX2=1 KX=1 C open(1,file='SPIDER_3.IN') C read(1,*) NX41,NX2,NPAR,STEP0,EPS,KX,IAType,Fl_thr,are_arc,dep_arc C C read(1,*) E_min, E_max, del_min, del_max, T_fix C C C* NXi - number of Xi data, C* NPAR - number of parameters, C* STEP - initial step, C* EPS - accuracy parameter, C* KX=1 or 2 for displaying the 1st or the 2nd argument, C !read(1,*) (PARGuess(I),I=1,NPAR) ! Initial parameters C close(1) NPAR=N_dat_in(1) IAType=N_dat_in(2) Npic=N_dat_in(3) !NpicY=N_dat_in(4) Nfreq=N_dat_in(4) N_case_f=N_dat_in(5) N_case_d=N_dat_in(6) STEP0 =r_inp(1) EPS =r_inp(2) Fl_thr =r_inp(3) are_arc=r_inp(4) dep_arc=r_inp(5) E_min =r_inp(6) ParGuess(:)=ParmsIn(:,1) PARMIN(:)=ParmsIn(:,2) PARMAX(:)=ParmsIn(:,3) E_max=ParGuess(6) T_fix=ParGuess(7)*1d6 ! If(npar.lt.(5.6)) PARGuess(6)=5. C write(*,'('' NX1 NX2 NPAR STEP EPS KX''/ C /3I5,1P,2E10.3,0P,I5)') NX41,NX2,NPAR,STEP,EPS,KX C write(*,'('' Initial parameters:'',1P,5E10.3)') C &(PARGuess(I),I=1,NPAR) C C PARMIN(1)=1.E-4 !Non-thermal density C PARMAX(1)=2.E3 C PARMIN(2)=1.E-2 !B-field C PARMAX(2)=3.E1 ! angle C PARMIN(3)= 22. !3. !0.07 ! 1.E-4 C PARMAX(3)= 87. !135 !177. !8. !1.E2 !-1 CCC** PARMIN(4)=4.E0 CCC** PARMAX(4)=8.E0 C CCC** PARMIN(5)=1.E-2 CCC** PARMAX(5)=1.E4 C PARMIN(4)=1e-2 !2.E0 Thermal density C PARMAX(4)=6e2 !88. !1.78E2 C ! PARMIN(5)=1.10E-1 !delta1 C ! PARMAX(5)=2.E1 C PARMIN(5)=del_min !1.10E-1 !delta1 C PARMAX(5)= del_max !2.E1 C PARMIN(6)=0.1 !0.003 !4.E0 E_max C PARMAX(6)= 10. ! 0.03 !20. !0.5 !176.E0 !2 C PARMIN(7)=1.5E-1 C PARMAX(7)=6.E0 !2 C PARMIN(8)=0.2E0 C PARMAX(8)=16.E0 !2 do 1 I=1,NPAR PSCALE(I)=ParGuess(I)/2. 1 continue 10 continue Cc Start of the fitting: select '0' for automatic fitting - standard mode. Cc Later change this to the default value Cc print*,'Input number of the manually varying parameter', Cc *' or 0 for automatic fitting:' Cc read*,KPAR Kpar=0 LMAX=1 C if (KPAR.ne.0) then C LMAX=0 C write(*,'('' Parameter ('',I2,''):'')') KPAR C read*,PAR(KPAR) C endif C Pause !Open(8,File='C:\GS_Modeling\Gelu_Blind_Test_1\DATA\Dim_DATA.dat') C Open(8,File='Dim_DATA.dat') C Read(8,*) NpicX,NpicY,Nfreq C !print*,'NpicX,NpicY,Nfreq=',NpicX,NpicY,Nfreq C !Pause C close(8) C C Allocate ( FreqData(1:Nfreq), C &Flcp(1:Npic,1:Nfreq),Frcp(1:Npic,1:Nfreq), C &ErrFlcp(1:Npic,1:Nfreq),ErrFrcp(1:Npic,1:Nfreq) C &,ArrNrl(1:Npic),ArrBgauss(1:Npic) C &,ArrEmax(1:Npic),ArrAngle(1:Npic) C &,ArrNth(1:Npic),ArrDelta1(1:Npic) C &,ErrNrl(1:Npic),ErrBgauss(1:Npic) C &,ErrEmax(1:Npic),ErrAngle(1:Npic) C &,ErrNth(1:Npic),ErrDelta1(1:Npic) C &,ErrData(1:Npic),RChi2(1:Npic) C & ,FFitR(1:Npic,1:Nfreq),FFitL(1:Npic,1:Nfreq) C &,ArrT(1:Npic),ErrT(1:Npic) ) Cc &,ArrParms(1:NpicX,1:NpicY,8),ErrParms(1:NpicX,1:NpicY,8)) !Assigning zeros to all arrays ArrNrl=0d0 ArrBgauss=0d0 ArrEmax=0d0 ArrAngle=0d0 ArrNth= 0d0 ArrDelta1=0d0 !If(Npar.lt.6) ArrDelta1=E_max ErrNrl= 0d0 ErrBgauss=0d0 ErrEmax=0d0 ErrAngle=0d0 ErrNth=0d0 ErrDelta1=0d0 ErrData=0d0 RChi2=0d0 FFitR=0d0 FFitL=0d0 ArrT=0d0 !If(Npar.lt.7) ArrT=T_fix ErrT=0d0 !End of Assigning zeros to all arrays !CALL AllDataSet(freqdata,Flcp,Frcp,ErrFlcp,ErrFrcp) !This reads arrays of the distribution function and !This reads arrays of the distribution function and !its parameters E,S,Mu,t freqdata=fr_in Flcp(:,:)=spec_in(:,:,1) Frcp(:,:)=spec_in(:,:,2) ErrFlcp(:,:)=spec_in(:,:,3) ErrFrcp(:,:)=spec_in(:,:,4) C1212 Format(31(1x,G12.5)) !,1x,G12.5,1x,G12.5,1x,G12.5,1x,G12.5) C open(1,file='TimeProMu.dat') c Do Imu=1,60 c Write(1,1212)xMu(Imu),(Fit1Data(14,IMu,JS1),JS1=31,60) c End Do c close(1) * open(1,file='SPIDER.IN') * read(1,*) NX41,NX2,NPAR,STEP0,EPS,KX * NXi - number of Xi data, * NPAR - number of parameters, * STEP - initial step, * EPS - accuracy parameter, * KX=1 or 2 for displaying the 1st or the 2nd argument, * read(1,*) (PAR(I),I=1,NPAR) ! Initial parameters * close(1) NYmin=1 !17 !20 !23 !28 NYmax=Npic !17 !20 !28 !29 NXmin=1 !10 !1 !NpicX-2 !16 !1 !NpicX-3 !1 !21 !15 !1 !6 !30 NXmax=Npic !NXmin+10 !NpicX !11 !6 !30 !60 Nx1=Nfreq !64 !30 ! !write(*,*) 'X2= ', X2dat C open(22,file='TErrFastElDens.dat') C open(3,file='TErrBgauss.dat') C open(4,file='TErrTheta.dat') C open(5,file='TErrThElDens.dat') C open(15,file='TErrDelta1.dat') C open(25,file='TErrEmax.dat') C open(27,file='XChi2.dat') C open(24,file='TErrT.dat') Par=ParGuess IparCorrX=Npic !100 !102 !18 IparCorrY=Npic !100 C Do JMu=NMumin, NMumax !25,NMumax ! NMumin, NMumax !36, NJMu !Nx=1,NxTot 1,NJS !32,32 ! C ! !$OMP PARALLEL DO Do Nx=NXmin, NXmax !,4 !12, NEmax !NEmin, NEmax !NJE,NJE !1, 29 !NJE,NJE 28,28 !Do Ny=NYmin, NYmax !,4 !21,NSmax ! NSmin, NSmax !Do Nx=NXmin, NXmax !,4 !12, NEmax !NEmin, NEmax !NJE,NJE !1, 29 !NJE,NJE 28,28 ! !If(Nx.lt.IparCorr.OR.Nx.eq.60.OR.Nx.eq.64) then If(Nx.le.IparCorrX) then Par=ParGuess DPar=ParGuess EndIf Step=Step0 * print*,Par,Step * Pause call SpecPixel(Nx,freqdata,Flcp,Frcp,ErrFlcp,ErrFrcp) !Call the time profile to fit with !the trapping function TestFlux=sum(Fundat) TestFluxR=sum(Fundat(:,1)) TestFluxL=sum(Fundat(:,2)) StocksV=TestFluxR-TestFluxL If(StocksV.lt.0) Par(3)=180.-Par(3) !Change trial viewing angle for L polarization !print*,'TestFlux=',TestFlux,'Nx #',Nx,'Ny #',Ny !If(TestFlux.le.22.OR.(TestFlux.ge.23.AND.TestFlux.le.32)) then ! 40 for unfolded model If(TestFlux.le.(Fl_thr)) then ! 40 for unfolded model !If(TestFlux.le.22) then ! 40 for unfolded model !print*,'TestFlux=',TestFlux !print*,'TestFlux=',TestFlux,'Nx #',Nx,'Ny #',Ny go to 2104 End If CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC print*,'TestFlux=',TestFlux,'Nx #',Nx,'Ny #',Ny !Pause Cc the fitting is being done by calling the SIMPLEX function ix1=0 do 4 L=0,LMAX Cc write(*,*)Fundat(1:10,1) if (L.eq.1) call SIMPLEX AvErr=Delta(ErrMx,IX1,IX2,RedChi2) X1=X1dat(IX1) !X2=X2dat(IX2) C if (NPAR.gt.0) then C C write(*,'('' Parameters:''/(1P,1X,10E11.4))') (PAR(I),I=1,NPAR) C write(*,'('' DeltaParam:''/(1P,10E11.4))') (DPAR(I),I=1,NPAR) C endif CCc write(*,*)Fundat(1:10,1) C C C write(*,'('' RMS error:'',1PE11.4)') AvErr C print*,' ' C Write(*,*)'Chi2=',RedChi2 C print*,' ' C write(*,'('' MAX error:'',1PE11.4,'' at X1 ='',E10.3,'', X2 ='', C *E10.3)') ErrMx,X1,X2 !Pause 4 continue C open(1,file='SPIDER.RES') C if (NPAR.gt.0) then C write(1,'(''Parameters:''/(1P,10E11.4))') (PAR(I),I=1,NPAR) C write(1,'(''DeltaParam:''/(1P,10E11.4))') (DPAR(I),I=1,NPAR) C endif C write(1,'(''RMS error:'',1PE11.4)') AvErr C write(1,'(''MAX error:'',1PE11.4,'' at X1 ='',E10.3,'', X2 ='', C *E10.3)') ErrMx,X1,X2 C close(1) C call DRAWING(KPAR) !open(9,file='C:\GS_Modeling\Gelu_Blind_Test_1\DATA\FitRes_2011.dat') 612 Format(1x,G12.5,6(1x,G12.5)) C open(9,file='FitRes_2018.dat') CCcc writing the original and fitted light curves for a specific set of parameters E,Mu,and S Do I=1,NX1 xxxMu=X1dat(I) Dum=0 FitProf=FitFun(xxxMu,dum) !/100. CCc Spectrum=Fit1Data(NE,I,JS) !FunDat(I,1) ! CCc ErrProf=Err1Data(NE,I,JS) CCc CCcc Write(9,12)xxxMu,TimeProf,ErrProf,FunDat(I,1),FitProf Cc Write(9,612)xxxMu,FunDat(I,2),ErrDat(I,2),dum , Cc &FunDat(I,1) ,ErrDat(I,1),FitProf C Write(9,612)xxxMu,FunDat(I,2)+FunDat(I,1),ErrDat(I,2)+ C &ErrDat(I,2),dum+FitProf Cc &,FunDat(I,2),ErrDat(I,2) !,dum FFitR(Nx,I)=FitProf FFitL(Nx,I)=dum C End Do C Close(9) Ccspi Pause ccspi w=w*1.1 ccspi End Do Ne=Nx Js=Ny ErrData(NX)=AvErr ArrNrl(NX)=1e7*Par(1) !Fast electron number density ArrBgauss(NX)=Par(2)*1e2 !Absolute value of magnetic field ArrEmax(NX)=Par(3) !Theta in Degrees ArrAngle(NX)=Par(4)*0.1e10 !Thermal density ArrNth(NX)=Par(5) !delta1 ArrDelta1(NX)=Par(6) !Emin ArrT(Nx)=Par(7) RChi2(NX)=RedChi2 If(Npar.lt.6) ArrDelta1(NX)=E_max If(Npar.lt.7) ArrT(NX)=T_fix !!!Use of AvErr is needed to form 'weights' of the data points for further fitting or analysis !ErrNrl(NX,NY)=1e7*sqrt((AvErr*Par(1))**2+DPAR(1)**2)/4. !*1e16 !ErrBgauss(NX,NY)=sqrt((AvErr*Par(2))**2+DPAR(2)**2)*1e2/4. !ErrEmax(NX,NY)=sqrt((AvErr*Par(3))**2+DPAR(3)**2)/4. !*7.27e7 !ErrAngle(NX,NY)=sqrt((AvErr*Par(4))**2+DPAR(4)**2)*0.1e10/4. !Thermal density !ErrNth(NX,NY)=sqrt((AvErr*Par(5))**2+DPAR(5)**2)/4. !*1e10 !ErrDelta1(NX,NY)=sqrt((AvErr*Par(6))**2+DPAR(6)**2)/4. !*1e2 !Use DPAR only below to estimate fit errors only CorrErr=1. ErrNrl(NX)=1e7*DPAR(1)/CorrErr !*1e16 ErrBgauss(NX)=DPAR(2)*1e2/CorrErr ErrEmax(NX)=DPAR(3)/CorrErr !*7.27e7 ErrAngle(NX)=DPAR(4)*0.1e10/CorrErr !Thermal density ErrNth(NX)=DPAR(5)/CorrErr !*1e10 ErrDelta1(NX)=DPAR(6)/CorrErr !*1e2 ErrT(Nx)=DPar(7) If(Npar.lt.6) ErrDelta1(NX)=0d0 If(Npar.lt.7) ErrT(NX)=0d0 Cc++++++++++++++below we consider the case of symmetric magnetic trap: integration time is reduced Cc by the factor of 2: JMu changes 1 to 30 rather than 1 to 60 2104 Continue Xcm=2.*(Nx-1) !Arcseconds along Y axes C If (ArrBgauss(Nx,Ny).gt.1e-2) then C C Write(22,*)Xcm,ArrNrl(Nx,Ny),ErrNrl(Nx,Ny) !Fast e density C Write(3,*)Xcm,ArrBgauss(Nx,Ny),ErrBgauss(Nx,Ny) !B field C Write(4,*)Xcm,ArrEmax(Nx,Ny),ErrEmax(Nx,Ny) !Theta (viewing angle) C C Write(5,*)Xcm,ArrAngle(Nx,Ny),ErrAngle(Nx,Ny) !Thermal density C Write(15,*)Xcm,ArrNth(Nx,Ny),ErrNth(Nx,Ny) !delta1 C Write(25,*)Xcm,ArrDelta1(Nx,Ny),ErrDelta1(Nx,Ny) !Emin C Write(27,*)Xcm,RChi2(Nx,Ny) C Write(24,*)Xcm,ArrT(Nx,Ny),ErrT(Nx,Ny) C C End If DPar(:)=Par(:) If(Nx.le.IparCorrX) DPar=ParGuess C Do i=1,Nx1 C write(*,*)NE !X1dat(i),x1 C end do !If(TestFlux.gt.(3)) pause End Do !Nx C End Do !Ny C ! !$OMP END PARALLEL DO C Close(22) C Close(3) C Close(4) !Close(5) c Close(15) C Close(25) C Close(27) C Close(24) c open(3,file='c:\GS_Modeling\EOVSA_Model_Fit_Test_1 c &\Some_Fit_Results\AlSp.dat') spec_out(:,:,1)=FFitL(:,:) spec_out(:,:,2)=FFitR(:,:) C open(3,file='AllFitSpecR.dat') C open(4,file='AllFitSpecL.dat') C Write(3,*)FFitR(:,:,:) C Write(4,*)FFitL(:,:,:) C Close(3) C Close(4) C C C open(2,file='ArrNrl.dat') C open(3,file='ArrBgauss.dat') C open(4,file='ArrTheta.dat') C open(5,file='ArrDensTherm.dat') C open(15,file='ArrDelta1.dat') C open(25,file='ArrEmin.dat') C open(24,file='ArrT.dat') C open(27,file='ArrChi2.dat') C C Write(2,*)ArrNrl(:,:) C Write(3,*)ArrBgauss(:,:) C Write(4,*)ArrEmax(:,:) C C Write(5,*)ArrAngle(:,:) C Write(15,*)ArrNth(:,:) C Write(25,*)ArrDelta1(:,:) C Write(24,*)ArrT(:,:) C Write(27,*)RChi2(:,:) C C C Close(2) C Close(3) C Close(4) C Close(5) C Close(15) C Close(25) C Close(24) C ArrParms(:,1)=ArrNrl(:) !Non-thermal density*1d7 ArrParms(:,2)=ArrBgauss(:) !B-field*1e2 ArrParms(:,3)=ArrEmax(:) !theta, degrees ArrParms(:,4)=ArrAngle(:) !Thermal density * 1d9 ArrParms(:,5)=ArrNth(:) !delta1 ArrParms(:,6)=ArrDelta1(:) !E_max, MeV ArrParms(:,7)=ArrT(:) !T *10^7 K ArrParms(:,8)=RChi2(:) !Fit Chi-2 ErrParms(:,1)=ErrNrl(:) !Non-thermal density*1d7 ErrParms(:,2)=ErrBgauss(:) !B-field*1e2 ErrParms(:,3)=ErrEmax(:) !theta, degrees ErrParms(:,4)=ErrAngle(:) !Thermal density * 1d9 ErrParms(:,5)=ErrNth(:) !delta1 ErrParms(:,6)=ErrDelta1(:) !E_max, MeV ErrParms(:,7)=ErrT(:) !T *10^7 K ErrParms(:,8)=ErrData(:) !Average fit error C open(2,file='ErrNrl.dat') C open(3,file='ErrBgauss.dat') C open(4,file='ErrTheta.dat') C open(5,file='ErrDensTherm.dat') C open(15,file='ErrDelta1.dat') C open(25,file='ErrEmin.dat') C open(24,file='ErrT.dat') C C C open(22,file='ErrData.dat') C C C Write(2,*)ErrNrl(:,:) C Write(3,*)ErrBgauss(:,:) C Write(4,*)ErrEmax(:,:) C C Write(5,*)ErrAngle(:,:) C Write(15,*)ErrNth(:,:) C Write(25,*)ErrDelta1(:,:) C Write(24,*)ErrT(:,:) C C Write(22,*)ErrData(:,:) C C C Close(22) C C C Close(2) C Close(3) C Close(4) C Close(5) C Close(15) C Close(25) C Close(27) C Close(24) !Npic=Nx*Ny * Do Nx=NXmin, NXmax !12, NEmax !NEmin, NEmax !NJE,NJE !1, 29 !NJE,NJE 28,28 ! * Xcm=7.327e7*(Nx-1) * Do Ny=NYmin, NYmax !21,NSmax ! NSmin, NSmax * * If (ArrBgauss(Nx,Ny).gt.0) then * * Write(2,*)Xcm,ArrNrl(Nx,Ny),ErrNrl(Nx,Ny) * Write(3,*)Xcm,ArrBgauss(Nx,Ny),ErrBgauss(Nx,Ny) * Write(4,*)Xcm,ArrDepth(Nx,Ny),ErrDepth(Nx,Ny) * * Write(5,*)Xcm,ArrAngle(Nx,Ny),ErrAngle(Nx,Ny) * Write(15,*)Xcm,ArrNth(Nx,Ny),ErrNth(Nx,Ny) * Write(25,*)Xcm,ArrDelta1(Nx,Ny),ErrDelta1(Nx,Ny) * End If * End Do * End Do Cc Do JS= 1,NJS !36,NJMu ! Cc Do NE=NEmin,NEmax !7,11 !NJE, NJE !1,29 !JE,JE !27,27 !1,1 !29 !Write(2,12)Energy(NE),TauR(NE,JS),ErrTau(NE,JS) !,Err(NE,JS)*TauR(NE,JS) !Alldata(Nx,Ny,20) ! !Ampl(Ltime) !Write(3,12)Energy(NE),Ampl(NE,JS),ErrAmpl(NE,JS) !,Err(NE,JS)*Ampl(NE,JS) !Write(2,12)Scm(JS),TauR(NE,JS) !Alldata(Nx,Ny,20) ! !Ampl(Ltime) !Write(3,12)Scm(JS),Ampl(NE,JS) Cc Write(2,12)xMu(JS),TauR(NE,JS),ErrTau(NE,JS) !Alldata(Nx,Ny,20) ! !Ampl(Ltime) Cc Write(3,12)xMu(JS),Ampl(NE,JS),ErrAmpl(NE,JS) Cc End Do Cc End Do 12 Format(1x,G12.5,1x,G12.5,1x,G12.5,1x,G12.5,1x,G12.5) !,1x,G12.5,1x,G12.5,1x,G12.5 &,1x,G12.5,1x,G12.5) !print*,'To continue? (Y/N)' !read(*,'(A1)') KEY !if (KEY.ne.'N'.and.KEY.ne.'n') goto 10 return !stop end function RESIDUAL(DAT1,ERR1,FIT1,DAT2,ERR2,FIT2) ! The norm which square is minimizing IMPLICIT REAL*8 (A-H,O-Z) common/Weight/ErrMean,FunMean,WeightMean common/Fit_case/N_case_d,N_case_f !RESIDUAL=WeightMean*(DAT-FIT)**2/(ERR+0.4e-1*DAT+1.)**2 If(N_case_d.eq.1) N_case_f=1 If(N_case_d.EQ.1) then If(N_case_f.EQ.1) then !To fit Stokes I, one has to use two lines below: RESIDUAL=WeightMean*(DAT1-FIT1+DAT2-FIT2)**2/((ERR1 & +ERR2)**2) EndIf EndIf If(N_case_d.EQ.2) then If(N_case_f.EQ.1) then !To fit Stokes I, one has to use the line below: RESIDUAL=WeightMean*(DAT1-FIT1+DAT2-FIT2)**2/(ERR1**2+ERR2**2) EndIf If(N_case_f.EQ.2) then RESIDUAL=WeightMean*((DAT1-FIT1)**2/(ERR1)**2 & +(DAT2-FIT2)**2/(ERR2)**2) !Now the polarized (R&L) spectra are fit. EndIf If(N_case_f.EQ.3) then RESIDUAL=WeightMean*( & (DAT1-FIT1+DAT2-FIT2)**2/(ERR1**2+ERR2**2) & +(DAT1-FIT1-DAT2+FIT2)**2/(ERR1**2+ERR2**2)) !Now the polarized spectra (I&V) are fit EndIf If(N_case_f.EQ.4) then RESIDUAL=WeightMean*(DAT1-FIT1+DAT2-FIT2)**2/((ERR1)**2 & +(ERR2)**2) &+((DAT1-DAT2)/(DAT1+DAT2+1d-5)-(FIT1-FIT2)/(FIT1+FIT2+1d-5))**2/ &(((ERR1+1d-5)/(DAT1+1d-5))**2 +((ERR2+1d-5)/(DAT2+1d-5))**2) c RESIDUAL=WeightMean*((DAT1-FIT1)**2/(ERR1)**2+ c &(DAT2-FIT2)**2/(ERR2)**2 c &+((DAT1-DAT2)/(DAT1+DAT2+0.1)-(FIT1-FIT2)/(FIT1+FIT2+0.1))**2/ c &(ERR2/DAT2)**2) !Now the polarized spectra (I&P) are fit EndIf EndIf return end !!!!!!!!!!!!!!! The program calls SPIDRSUB !!!!!!!!!!!!!!!!!!!!!!!!!!!!!