dynapav.cpp

/*Copyright Cheng Sheng 2007*/
#include "arprec/mp_real.h"
#include "arprec/mp_complex.h"
#include "arprec/mp_int.h"
#include <fstream>
#include <iostream>
#include <time.h>
//#include <sys/time.h>

using namespace std;
using std::cin;
using std::cout;
using std::endl;
	
#define   EPS   1E-50  
#define   TRUE     1  
#define   FALSE   0 


int GSolve(mp_complex **a,int n,mp_complex *x) 
{ 
	int i,j,k,maxrow; 
 	mp_complex tmp,*a2;
	a2=new mp_complex[n];	
		
	for (i=0;i<n;i++) { 
 		/* Find the row with the largest first value */ 
 		maxrow = i; 
 		for (j=i;j<n;j++) { 
 			if (abs(a[j][i]) > abs(a[maxrow][i])) 
  			maxrow = j; 
		} 
 
 		/* Swap the maxrow and ith row */ 
 		for (k=i;k<n;k++) { 
  			tmp = a[i][k]; 
  			a[i][k] = a[maxrow][k]; 
  			a[maxrow][k] = tmp; 
 		} 
		tmp=x[i];
		x[i]=x[maxrow];
		x[maxrow]=tmp;
		 
 		/* Singular matrix? */ 
 		if (abs(a[i][i]) < EPS) 
  		return(FALSE); 
 
 		/* Eliminate the ith element of the jth row */ 
		tmp=a[i][i];
		for (k=i+1;k<n;k++) {a2[k]=a[k][i];}
		x[i]/=tmp;	
 		for (j=i;j<n;j++) 
		{
			a[i][j]/=tmp;					
  			for (k=i+1;k<n;k++) 
			{ 
  				a[k][j] -= a2[k] * a[i][j] ; 
  			} 
 		} 
		for (k=i+1;k<n;k++) 
		{ 
			x[k] -= a2[k]*x[i];
  		} 
 	} 
 
 	/* Do the back substitution */ 
	for (i=n-1;i>=0;i--) { 
 		tmp = mp_complex(0.0,0.0); 
 		for (j=i+1;j<n;j++) tmp += a[i][j] * x[j]; //modified
 		x[i] -= tmp ; 
 	} 
 
 	return(TRUE); 
}
 

mp_real bessj0(mp_real x)
	//Returns the Bessel function J0(x) for any real x.
{
	mp_real ax,z;
	mp_real xx,y,ans,ans1,ans2; //Accumulate polynomials in mp_real precision.
	if ((ax=abs(x)) < 8.0) { //Direct rational function
		y=x*x;
		ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
				+y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
		ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
				+y*(59272.64853+y*(267.8532712+y*1.0))));
		ans=ans1/ans2;
	} else { //Fitting function (6.5.9).
	z=8.0/ax;
	y=z*z;
	xx=ax-0.785398164;
	ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
			+y*(-0.2073370639e-5+y*0.2093887211e-6)));
	ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			+y*(-0.6911147651e-5+y*(0.7621095161e-6
			-y*0.934935152e-7)));
	ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
	}
	return ans;
}

int main(int argc, char *argv[])
{
	mp::mp_init(50);
	mp::mpsetprecwords(50);
	
	mp_complex I;
	I=mp_complex(0.0,1.0);
	mp_real Pi=3.1415926535897932384626433832795;
	const int M=300;
	const int N=400;
	
	//struct timeval tpstart,tpend; 
	//float timeuse; 
	clock_t start,stop;
	
	double P0,R,T0,Tt,t0,a;
	double *points;
	double *E,*nu,*rho,*h,*t;
	int *friction,*layerindex,nr_layerindex;
	mp_real *mu,*lambda,*cp,*cs,*G;
	mp_complex **Ksys,***w,*P;
	int nr_sys;
	int sys_nr;
	int layers,nr_points;
	int i,j,layers_nr,points_nr,n,m;
	mp_complex K[4][4],K0[2][2];
	mp_complex s,ts,sin_s,cos_s,sin_ts,cos_ts,DD;
	mp_complex s0,ts0;
	mp_complex Ft[N],Fn[N];
	mp_real omega[N];
	mp_complex W[N],T[M],Wt[N];
	mp_real **Result,**Figure;
	mp_real Min,time_Min;
	mp_real r;
	//double temp;

	
	const mp_real alpha[]={2.4048255576957727686216318793264546431242449091460,
5.5200781102863106495966041128130274252218654787829,
8.6537279129110122169541987126609466855657952312754,
11.791534439014281613743044911925458922022924699695,
14.930917708487785947762593997388682207915850115633,
18.071063967910922543147882975618176560248986747001,
21.211636629879258959078393350526306836181808975976,
24.352471530749302737057944763178907184569372675149,
27.493479132040254795877288234607414546529568860550,
30.634606468431975117549578926854232737273571629178,
33.775820213573568684238546346714728023932705742567,
36.917098353664043979769493063272952636485702585576,
40.058425764628239294799307373994472910347961951990,
43.199791713176730357524072728743422170876097920882,
46.341188371661814018685788879112849174650148505260,
49.482609897397817173602761533178272269702237753039,
52.624051841114996029251285380391573300120332932687,
55.765510755019979311683492773461830631375339755205,
58.906983926080942132834406634615685585658088855756,
62.048469190227169882852500264650952323815383766251,
65.189964800206860440636033742512316225656567704777,
68.331469329856798270992303839984378036847621222727,
71.472981603593732825063073856129700553371893825789,
74.614500643701837883820540469335554203020853784373,
77.756025630388055037739371891233693265083562216138,
80.897555871137627863772143490873052856874725969662,
84.039090776938190157879638347997955127712589417406,
87.180629843641153651261805069048608848672129946456,
90.322172637210480055717766777627908914486121572578,
93.463718781944774171190591543980910757323940062918,
96.605267950996268778121617323927896879425560391421,
99.746819858680596470279979000134577801039929537751,
102.88837425419479459642003427256457628927470004696,
106.02993091645161551017691719188244758163597322099,
109.17148964980538355206597701274965661829727689521,
112.31305028049490962749450612217833138594062200411,
115.45461265366693962811775669402959659734658011560,
118.59617663087253171562938447524390966025681886518,
121.73774208795096296523436348334089451029740458848,
124.87930891323294604525912836690414362854635369113,
128.02087700600832407976355939229794390008143742793,
131.16244627521391460789611499152693779453700440921,
134.30401663830546609935289683410038249518973723559,
137.44558802028427778778266148333216574501469751324,
140.58716035285429654848890334006056965682193987633,
143.72873357368973253395071800477048555830997691241,
146.87030762579664959413270485825331810328877264830,
150.01188245695475749088047297407295686286872420248,
153.15345801922789248759162561015357879452145238777,
156.29503426853352381954949527313604669838971432159,
159.43661116426314632349103791000180161279196617144,
162.57818866894667751905980896294183321651444715279,
165.71976674795502086669042673820362638640107738548,
168.86134536923582568745633745233179715962999675443,
172.00292450307820021540043150732572603521334713255,
175.14450412190274306537279785910344570599465017568,
178.28608420007377068148397419125075748636903921768,
181.42766471373105079421219958705073504225370515147,
184.56924564063871814111541224791702088292439432262,
187.71082696004935978007526572080084190428715565842,
190.85240865258152232178079894959733273054839191129,
193.99399070010911978994039672748625030235183028735,
197.13557308566141473621200268072495795270830997233,
200.27715579333241178336208747126910756890452518549,
203.41873880819864617124882021531933310925718874770,
206.56032211624447365545724492840174638521043239287,
209.70190570429407519748244666580127306056180741375,
212.84348955994948275074556567746795366135512841016,
215.98507367153401315699563788314249439745514051660,
219.12665802804056746518913097823472201395539836269,
222.26824261908431434128452437947225138877036639484,
225.40982743485932989851318199627000062429318381297,
228.55141246609881330119784069022879189489298203722,
231.69299770403853878097379619250103752347952276116,
234.83458314038324101980802326656358016257541394657,
237.97616876727566285555756332328575945849495437339,
241.11775457726802251499697108991101422808443230154,
244.25934056329568255880657617786596326564904335797,
247.40092671865282484808794068946977122248473266764,
250.54251303696995547044636528300125132156850424656,
253.68409951219308100469280850438583235415692045925,
256.82568613856441302431773328636516484357503628270,
259.96727291060447157306477756919990550063230282933,
263.10885982309547069270400257633712703627333307222,
266.25044687106588011884095423400708460818973985461,
269.39203404977606713817339613835350035429537058004,
272.53362135470493145351434896934894548640986317658,
275.67520878153745384789779496894639374397041426683,
278.81679632615308657846177386447988752366426589122,
281.95838398461491985431670967679016149868657908406,
285.09997175315956453913058774304653504916788468579,
288.24155962818769643811739502072485249215625997429,
291.38314760625521224168905006922144698187326414241,
294.52473568406495145823317711091537938571856031660,
297.66632385845894252404919067485237989556409082931,
300.80791212641113477166896885131883028206660019937,
303.94950048502058110605777975398170928572619805691,
307.09108893150503911477379700032604989002392778639,
310.23267746319496095259609304361402241300011466251,
313.37426607752784471969024510190451272294145129640,
316.51585477204292221847752267425878978745631761552,
319.65744354437615994888734721776411877350897639715,
322.79903239225555200324754484869982270800535646078,
325.94062131349668516737048895575563303978138729129,
329.08221030599855803834301990803369632112684181664,
332.22379936773963734548991829833928045685758204914,
335.36538849677413592094340034522326340219801714048,
338.50697769122849792097948562312190517150914740745,
341.64856694929807795845225875691782244429713122518,
344.79015626924400177897536108638471039854560841911,
347.93174564939019700680817750011965491366832819223,
351.07333508812058330777196349075980345298650288661,
354.21492458387641207232872856145917956570114621807,
357.35651413515374641796625372799147456998186727906,
360.49810374050107295145694339569459436282704668128,
363.63969339851703732311190421565425400750523572844,
366.78128310784829615110604419024673014152295828433,
369.92287286718747839817743556215557503805765047378,
373.06446267527124974901834601149118125410717537342,
376.20605253087847396766713604931075750630985015300,
379.34764243282846561307989991617750014290528180126,
382.48923237997932886044652575980009292193413142964,
385.63082237122637751811959420587623861286716902783,
388.77241240550063164743456431331437635091756329611,
391.91400248176738648721140418522647947445449147027,
395.05559259902484965816206392402355583213515887812,
398.19718275630284287644780017478839413178146966061,
401.33877295266156464175501052544932939332792816635,
404.48036318719041058487775971257864770540842297917,
407.62195345900684836418159264618835495296786475127,
410.76354376725534419064153695411882161002251629783,
413.90513411110633823846514189925424338540885461941,
417.04672448975526636360825131453658214832756960886,
420.18831490242162570666570138589883671758164544249,
423.32990534834808190049664257124831984719144308750,
426.47149582679961573728320348173278102963861577882,
429.61308633706270727521922826791818567145736886080,
432.75467687844455548232461790851459773924346559607,
435.89626745027233162457130126955958910124500189470,
439.03785805189246470813345549591331034811537230834,
442.17944868266995738163914409412261753442470032345,
445.32103934198773079426595862781050020361271472300,
448.46263002924599698981678026027863169137514750419,
451.60422074386165749592797226987966919736595259317,
454.74581148526772684166574391328349769972831399382,
457.88740225291277980629414575823244661255887787683,
461.02899304626042126726195909104513305583681698566,
464.17058386478877757674420737417214615408590653207,
467.31217470799000845365438654061521888074725472727,
470.45376557536983843216339407593589020618631480255,
473.59535646644710695865007811949473283691148636004,
476.73694738075333627687923248125502110594173451322,
479.87853831783231628625330535432252806029782027622,
483.02012927723970560039750756503156759412045812201,
486.16172025854264807328480056415728702237785378629,
489.30331126131940409774578087142879968206117971777,
492.44490228515899601668600916307383250631971670863,
495.58649332966086702078684377095595744370460885680,
498.72808439443455293802282359755481233794691210833,
501.86967547909936635010781928003145424114231174807,
505.01126658328409249909412547840127937840373979652,
508.15285770662669647389647235678633486200631664888,
511.29444884877404119159269543380744271792636309104,
514.43604000938161571205416906271737084050057567649,
517.57763118811327344686574343388851513681302921053,
520.71922238464097984468493957250457077617750040374,
523.86081359864456915523651010797574829833418429762,
527.00240482981150989310934034999316198640516551201,
530.14399607783667864048177529477162368310469978577,
533.28558734242214184490841475993210113011513239704,
536.42717862327694528441198256460414391752287023684,
539.56876992011691088739023694449615696264330410044,
542.71036123266444060931891545413996367238316697608,
545.85195256064832708195317202467976011686329355653,
548.99354390380357076374475894473724771606537575839,
552.13513526187120333254054012989642625212954564934,
555.27672663459811707334749969882382055476950090801,
558.41831802173690002507560380059492386191408999214,
561.55990942304567666073587936128739264700159886402,
564.70150083828795388560805507539313052414182784108,
567.84309226723247214742933517053045488086917612156,
570.98468370965306146172084152242081923168225129358,
574.12627516532850216398680588732201024690060944765,
577.26786663404239020871801572241781373437613823876,
580.40945811558300684292816220179805124896089282979,
583.55104960974319248937110251938625125888228770580,
586.69264111632022468164885924334740998107588049493,
589.83423263511569990014348104065821607502341793776,
592.97582416593541916410861541134146760505838932881,
596.11741570858927724135568942285772574535369462975,
599.25900726289115534278087885822423065823842490994,
602.40059882865881717451757647369716173262997235163,
605.54219040571380822577899659058820133166004974176,
608.68378199388135817549022304874291832201519685134,
611.82537359299028630561101379443901602207871938831,
614.96696520287290981363190358388155539684522080793,
618.10855682336495492109781748805940371108754683461,
621.25014845430547067918611895253088132668424172881,
624.39174009553674537634977748249752793264031256969,
627.53333174690422545684061172593907434186147303330,
630.67492340825643686256128465378469347979066410603,
633.81651507944490871416635205327631615129721209648,
636.95810676032409925065018840599065511835223051610,
640.09969845075132394983059855271991823302076375907,
643.24129015058668575516852617628727056531941471767,
646.38288185969300733726326526561872340765273073918,
649.52447357793576532113538215483132746372804566556,
652.66606530518302641306223766165446969375918019467,
655.80765704130538536326931541032135950227173377650,
658.94924878617590470320996598476564199268423809420,
662.09084053967005619849183490937272264591549205366,
665.23243230166566396073505493566666097532754010062,
668.37402407204284916377989575284830930208258827704,
671.51561585068397631170438397707011399158464033553,
674.65720763747360100806961431793381042414914541598,
677.79879943229841917768603720513173936814453383162,
680.94039123504721769399169554448570251131498017542,
684.08198304561082636685677004607813430820176707298,
687.22357486388207124728127527487082437940553561409,
690.36516668975572920703755374559170906718330577720,
693.50675852312848375282940877245974619327878362425,
696.64835036389888203599821206874682036547252363849,
699.78994221196729302020588906646936894426512310532,
702.93153406723586677086795722137242798630377303083,
706.07312592960849483139927491517897507489987443507,
709.21471779899077165257323365335972106224255087808,
712.35630967528995704248406118384546628800349939712,
715.49790155841493960574385553550421199901900465652,
718.63949344827620114164299364808257231245255567147,
721.78108534478578197205660875775941958983076653563,
724.92267724785724717089277031508852842402842381378,
728.06426915740565366785160285286169905086710086658,
731.20586107334751820020053596281412287574231923887,
734.34745299560078608717079692297599283819805221297,
737.48904492408480080244567561650306921863058641576,
740.63063685872027432104347164772008636937921512009,
743.77222879942925821769877144703367962338957202747,
746.91382074613511549461612955573816737193725315582,
750.05541269876249311721161276923127401351103198064,
753.19700465723729523717121971875926714363921888724,
756.33859662148665708284206776812022820313962221407,
759.48018859143891949763354717518139657673161011568,
762.62178056702360410774243270291359828247714816789,
765.76337254817138910112922110851456701783863558629,
768.90496453481408560026368983365786225797700371610,
772.04655652688461461172676540689546675281877979567,
775.18814852431698453630412628613590595776596052412,
778.32974052704626922373538192410432338131636445184,
781.47133253500858655679196849358713469387086297575,
784.61292454814107754984784745460162557633264539859,
787.75451656638188594758041789444041624983035009115,
790.89610858967013830989545725599001634429608749453,
794.03770061794592456961005716042835852712351165109,
797.17929265115027904985206197696031710371150847612,
800.32088468922516192854406039023749452380973693144,
803.46247673211344113773511090463776751235160216110,
806.60406877975887468592466231811865797793702242241,
809.74566083210609339189109705643870747443537109162,
812.88725288910058401889249339669867289141456870658,
816.02884495068867279845006475426874564401611138008,
819.17043701681750933325576230018653893616454108010,
822.31202908743505086906517332497348208299235313449,
825.45362116249004692574554515023516308635893311465,
828.59521324193202427794692812599850326215120414860,
831.73680532571127227615245920686335592604795708044,
834.87839741377882849914208117262324290918892701799,
838.01998950608646472917287742034145102884065151565,
841.16158160258667324143904905351396418057855116207,
844.30317370323265339962570602650129387013175684305,
847.44476580797829854961340995193977284852795949414,
850.58635791677818320362510232269704930010330729614,
853.72795002958755050733397528845813559366393479725,
856.86954214636229998267027773468910530750126907997,
860.01113426705897553927727076209053543242396900465,
863.15272639163475374777181633050548205160761207675,
866.29431852004743236816365291113601993984765634064,
869.43591065225541912697952456332466537906198445122,
872.57750278821772073682421740838200406901909712948,
875.71909492789393215229044332984886523515789310101,
878.86068707124422605630360942311709248781193957701,
882.00227921822934257115602937784677322389631815326,
885.14387136881057918864826768686029847494567944825,
888.28546352294978091391324971185555862851526066178,
891.42705568060933061765170319588431440533582751890,
894.56864784175213959165559572245476488882297579634,
897.71024000634163830263966703641632341511993212069,
900.85183217434176733954008772866153285592203021876,
903.99342434571696854957386299304425833281971236594,
907.13501652043217635848299247781101067472121252359,
910.27660869845280927051373945028182444052197032526,
913.41820087974476154380379387025834805999846877307,
916.55979306427439503696876857047395799374715612714,
919.70138525200853122279447458044156108151535580280,
922.84297744291444336505290488147252661955137501360,
925.98456963695984885456793508935087242820328866030,
929.12616183411290170076153982484139014632892398917,
932.26775403434218517501293568214982585662427378015,
935.40934623761670460226160247841365954392984174051,
938.55093844390588029738070667003777665874576264475,
941.69253065317954064294015348826338411491881683707};
	const mp_real Fm[]={0.00013357022854861339330752669743571005518850259255871,
0.00031089152528146017490991380138042664091171172161748,
0.00048839261427075016157294649209041388312424664371296,
0.00066578637866017244950715655386718541803939296450454,
0.00084300732230676509201898813381716214866545884792704,
0.0010200050247869966781632574884056948335359816472863,
0.0011967314896310576086375295962401168027429933465307,
0.0013731393621661448293173505059550766340182422166646,
0.0015491815549091772987599340971289477145239010621464,
0.0017248111521082802562220910897309482362417670123520,
0.0018999813866701549326975937849561810328738187576644,
0.0020746456397109410197669087619280825819369576537270,
0.0022487574482186954567883690535870837842264937459968,
0.0024222705159395103935468630371744181922774338325668,
0.0025951387256436060842838533107424210075818100927336,
0.0027673161520087485905743428950077335165685250898290,
0.0029387570747798365191731466451009591755644435870592,
0.0031094159920410772084995861791249514659168146382624,
0.0032792476335170364121310679482974337959085992312156,
0.0034482069738568317699360510426690573643382859156810,
0.0036162492458747090628769498141643757617134248788402,
0.0037833299537301276515089782535133924165291210232394,
0.0039494048860358372064168939527022960101152444818794,
0.0041144301288854325045466395168652762105235480273782,
0.0042783620787936151187936102335119049927259845647130,
0.0044411574555434335551895095713182237026460604930594,
0.0046027733149354179392485338732136768375045771864016,
0.0047631670614339375220162908459448495099777240475114,
0.0049222964607063828209598431437428575301462575761004,
0.0050801196520509634512913575744140380564285818743974,
0.0052365951607090498216442225705639765942881545202202,
0.0053916819100580913246280686962649195384759492874866,
0.0055453392336812274390322926040023698491256681629114,
0.0056975268873097787945691938693316867921251537252060,
0.0058482050606348675216316089859411784052146626109602,
0.0059973343889844732138942031692995523106681699123104,
0.0061448759648622845993229284523289913837619804565280,
0.0062907913493447588904966988938234856921161040317798,
0.0064350425833328516697170721457788433370851722672332,
0.0065775921986549306680351044723475889723214023651500,
0.0067184032290174373355734238351699826812090115632652,
0.0068574392207999109578992599379599369630170910592078,
0.0069946642436910414449593943511973902150041028203146,
0.0071300429011624689389435368026431630383421120654252,
0.0072635403407771011462777291821302253275607581447890,
0.0073951222643287728575040414767222002161879250321532,
0.0075247549378101265164023201249902759091164140871408,
0.0076524052012056479603123485961974216268684369828660,
0.0077780404781068475900417203700407508855056715604208,
0.0079016287851466342446642739874685867447332975633026,
0.0080231387412499869525639898921195900832395639657036,
0.0081425395766980884993892454782806781468704304870916,
0.0082598011420031443868730255100707864617207401973042,
0.0083748939165911712418589157722665062437202177162534,
0.0084877890172901000584428060846262760765977918817030,
0.0085984582066206018024022370726601209759442358939222,
0.0087068739008871058592823463879385782657377983138608,
0.0088130091780665455478622471253578948605318452748736,
0.0089168377854924294306245810933102608039206231165258,
0.0090183341473319024129835132540817479688296054882324,
0.0091174733718535266134933812764407910958716144423892,
0.0092142312584835786876574357876959022794272245433972,
0.0093085843046487276774403250506623252528078637508570,
0.0094005097124030255159290954474359225671967205491322,
0.0094899853948372110202106665608889187700190948533712,
0.0095769899822683975335537616188332705714647539783696,
0.0096615028282082843082329025978950141073422912284558,
0.0097435040151081022303859797316521048976433271129000,
0.0098229743598785755554920457649849631098513293142758,
0.0098998954191832529245093671821462336699056453843916,
0.0099742494945036330433380802962756221523532899769060,
0.010046019636974583008789485618221780313543748017387,
0.010115189651988620329200052734841066246342722332980,
0.010181744103567703193601961690771332585464470153532,
0.010245668318501247466178287043433421479652087801002,
0.010306948390249163198670476289795668803041666547707,
0.010365571182608778138415264326583918319687275969129,
0.010421524333144590739589867309667591944399965380324,
0.010474796256379870535745430254729372268726493371734,
0.010525376146749199378410568053858496678088093240804,
0.010573253981311122964953792654588222892644471144010,
0.010618420522220158244421883471963857268035876256476,
0.010660867318957478678057639710564142189880616873539,
0.010700586710319675916910966675325328043669908078519,
0.010737571826165073217595673809110899812438999082352,
0.010771816588917142823960749362140359430362944529342,
0.010803315714824656572344002192712129086248319012198,
0.010832064714978276106225752821683197869546990169589,
0.010858059896083366287540184996227974303966099322690,
0.010881298360988892641651418806571224708876453405671,
0.010901778008972340946067977104654164264531496593836,
0.010919497535780674344357311115162378954386661702606,
0.010934456433427420611440539468872016916157794961636,
0.010946654989746059389518823724069989848825295050262,
0.010956094287699956330350583108779046516642758372798,
0.010962776204449168094536465234049713595300477494972,
0.010966703410174519046987983464159364679273382929401,
0.010967879366659427225013480443624964146750717216663,
0.010966308325630033716663369671523621427648655874929,
0.010961995326854265947409705078490175346105320732283,
0.010954946196000541508237900932307996271177635722958,
0.010945167542256895043224450866360903001530105932614,
0.010932666755711386325174140378624342946230997299250,
0.010917452004494722959498703238653696034306929389946,
0.010899532231686106144944163042935188709181551712284,
0.010878917151983382560835052533842905335620319815694,
0.010855617248138659720137872096754783530714588794591,
0.010829643767160616001917710212666724607383334169831,
0.010801008716284810031869538871783374793265593334962,
0.010769724858713367091890209555590362990777878151495,
0.010735805709125492785609507870237101857947812061303,
0.010699265528960336243067373844097347231294593919956,
0.010660119321473796691123496913246855174867545339543,
0.010618382826570938223701810782078582266989422847443,
0.010574072515415748054773254591035896376271436128746,
0.010527205584820043404420757936924111644972659704392,
0.010477799951413401431961287103427094313661330607707,
0.010425874245596055267674145212882017621675470252610,
0.010371447805276767184165650136350555755543217539218,
0.010314540669397757267968028936588279955485643947964,
0.010255173571248832580029303331914630895356101265814,
0.010193367931572927708936808539767365497629323841298,
0.010129145851465332801903741250061509060560111311030,
0.010062530105068949584854854092669417230890698691691,
0.0099935441320679795337451207376575120078386195529272,
0.0099222120299825112141635590783370577845342225280906,
0.0098485585462665358452083789232574829833877605382500,
0.0097726090702119813467360063837326807372984087957660,
0.0096943896246614154768304058259042807767291235556014,
0.0096139268575321281394401927083664894161057938855596,
0.0095312480331543615216099366776601319885055042191612,
0.0094463810234265143869064995518310391211009134245710,
0.0093593542987902035881320859673742347709906859858848,
0.0092701969190281216501528630556118039077354128060824,
0.0091789385238876840949004583426719036590175330247040,
0.0090856093235335140178890291416552263747263225548954,
0.0089902400888318642618249469756524117722073660895060,
0.0088928621414701293512934105361695502499064152812740,
0.0087935073439146501366477920752193772844024015600478,
0.0086922080892100638290072018870499603958966561789166,
0.0085889972906235007759374583196466506696510502142726,
0.0084839083711369769135551653235651236892665077917188,
0.0083769752527913773204166372971105414687325763086274,
0.0082682323458854716769584550728861384151612911317140,
0.0081577145380334466871374702881712488121672384285722,
0.0080454571830844836323401863963499058234005753534284,
0.0079314960899079511880363082685158601439472098629222,
0.0078158675110478244278614864020792018268511999040934,
0.0076986081312499805550376754588846179379573282227456,
0.0075797550558660603248731416115753274594719020101228,
0.0074593457991376213425182934608532579164752107136324,
0.0073374182723643454255755645900312363994448519131002,
0.0072140107719600970003597203930928345966875460782894,
0.0070891619674006630427717409734772625588380998500206,
0.0069629108890670373694851938188996070492516524574480,
0.0068352969159881431224602116651768653587942516118490,
0.0067063597634869170601225657299638465653927002486074,
0.0065761394707337077627169113298760428243031931601244,
0.0064446763882109670686273173229101621280190203856810,
0.0063120111650932399745408560018578491754799627023410,
0.0061781847365464828473230673281501503455305460398080,
0.0060432383109507631019172893655297480371300367666600,
0.0059072133570504154904432883248568882187113936604894,
0.0057701515910357508163569879086786452628459553870744,
0.0056320949635604322278796284538222115276913127860194,
0.0054930856466986522511850963106915000018704273219010,
0.0053531660208462603907603936231385821253877177446410,
0.0052123786615700064470780407436533844706128578957246,
0.0050707663264090786758336823422650614717485734141664,
0.0049283719416331285345429736735435154294811556796000,
0.0047852385889609850277383975178163815141366529160246,
0.0046414094922442715682832390009317349530270992147238,
0.0044969280041201468167964230117594520067003862630300,
0.0043518375926373981416745842637555949509122418684402,
0.0042061818278601221569698456823507937969678556339536,
0.0040600043684532312431439734282298683676899430173156,
0.0039133489482540280356301778958070633013669024497958,
0.0037662593628340915777989679229094032048556080505926,
0.0036187794560557191783980598645032874575495899249306,
0.0034709531066271669893199409213418911977942035409554,
0.0033228242146609299285931490118655101515153979506904,
0.0031744366882392978171903763636567733330662615847116,
0.0030258344299914194784400949165498500758954593186518,
0.0028770613236861000678059945269569888793026760376516,
0.0027281612208445490612921652085377177437023196291410,
0.0025791779273772871359172276366261310938756881650924,
0.0024301551902494096291954514435823852717848830257052,
0.0022811366841783923704262302380382725287247973697590,
0.0021321659983686124393235523351451976814705143007058,
0.0019832866232867418320506456419772927952737122890684,
0.0018345419374821561064342179814197424303463262662536,
0.0016859751944564828428242787482821447095806893622909,
0.0015376295095863962009781530070129633817582442251335,
0.0013895478471037439831496378495571805056697707844868,
0.0012417730071370724363525219329826497758677535031048,
0.0010943476128185915500631886391870513880600942798624,
0.00094731409746059983737349555682970092012146971296972,
0.00080071469180536253616599267234807169502737998983344,
0.00065459141135241084103900784520861630647957062025402,
0.00050898604376720218565312118975515820112630712697468,
0.00036394013637505274850755239687383489455107802716666,
0.00021949498374422326289671737022000127107057900346432,
0.000075691615362007884359183560488617176021551719006212,
-0.000067429216592356682870812552911353469244937901847280,
-0.00020982704937417836385300912807726254399583532021292,
-0.00035146172169544472531648078442272927921060035789554,
-0.00049229338563889350841949757028256840933415374345440,
-0.00063228251846034370112212501629294832162987978040766,
-0.00077138993427565066704600963496842354339320155018648,
-0.00090957679562868849662420930722945131966416157562822,
-0.0010468046249368035330484957285192135527082912321668,
-0.0011830353158102249361292669265602070636186943971367,
-0.0013182311442419611674588184451340524176679866147533,
-0.0014523547796647553956195870749814269464951942185210,
-0.0015853692958717180156600032288568704429694198963978,
-0.0017172381817973007373792807377244183098410666890060,
-0.0018479253521553240064841382701274749161021324855927,
-0.0019773951579308178654334648435459290396436357848318,
-0.0021056123967224857204681698257155326625665152023710,
-0.0022325423229326508413037358778301879511281414010510,
-0.0023581506578015967632927586359158796435416399211188,
-0.0024824035992832650712797607350812164043535040072474,
-0.0026052678317593273022976948486318712470381074431610,
-0.0027267105355887018928187391093900251933706573022826,
-0.0028466993964896421972982488393032723600050738441020,
-0.0029652026147515775997762576334733021632724433951156,
-0.0030821889142739466105773105224332242489672965891990,
-0.0031976275514293185666499504455229967665832687016136,
-0.0033114883237481591175132421549317294871545365703056,
-0.0034237415784226540595652466093461993970244796878124,
-0.0035343582206270662598344262432561724303840541141922,
-0.0036433097216521613660444182024887062610056678931396,
-0.0037505681268512997127948175799675441506226770384598,
-0.0038561060633958542831762407747770027039774962799282,
-0.0039598967478376777504458661216705905071281134295914,
-0.0040619139934764054844739107529124455886034534822838,
-0.0041621322175294459412983263707302235339048888022830,
-0.0042605264481025750398499014218658805631632869087468,
-0.0043570723309591169460849224367141385721018824287214,
-0.0044517461360857601095432454760183256240403574326170,
-0.0045445247640531244087027205983708324429559325688686,
-0.0046353857521692628372114026313433749138799495820696,
-0.0047243072804243492807576314300167656115532022396658,
-0.0048112681772248725714288596915810659025100988896082,
-0.0048962479249157261401977763332695267038349883573458,
-0.0049792266650886521957918633926030457019973072525724,
-0.0050601852036755694166390119326674434860308270438870,
-0.0051391050158253836286897802825442902518637951108994,
-0.0052159682505629518324201789670216251660982427016346,
-0.0052907577352289412138204284772124439187233386997890,
-0.0053634569796993964031647027281752108651922318986656,
-0.0054340501803839002082188296082349333301902650835814,
-0.0055025222240012853215641484935994116218864832170942,
-0.0055688586911319270610934696334142340859833028086464,
-0.0056330458595457200245848613471492274084075302502446,
-0.0056950707073049145996225284483195644544864434932676,
-0.0057549209156410625449872101941458852271030117599140,
-0.0058125848716053943248993081523795416638785999250580,
-0.0058680516704920245090344465913211206224714374361640,
-0.0059213111180334553248695640660381388109299329215552,
-0.0059723537323679223404502377274934812168635171022476,
-0.0060211707457782002408631800888686807160554183514740,
-0.0060677541062015607163002975298164552764527329056936,
-0.0061120964785106485793511085751273080268523150688876,
-0.0061541912455651163497956321789207023514518508741202,
-0.0061940325090339316624332342724411773483286091058112,
-0.0062316150899883459431317617124983258546893427568590,
-0.0062669345292655868360952365933661759406249416396172,
-0.0062999870876034108271373115960263177026948098502192,
-0.0063307697455457263693597048521867281927556977656158,
-0.0063592802031195715549642405965917723192027717230922,
-0.0063855168792838039659223630346326488985859197536190,
-0.0064094789111499337528975806684017758991680249009262,
-0.0064311661529756042122447285916186556387188073175824,
-0.0064505791749312971312535532370513088689005629635056,
-0.0064677192616409129283069772656252258724893477675928,
-0.0064825884104969481036240615434638722713501475607632,
-0.0064951893297510647141930192144280730880249318825852,
-0.0065055254363809184699185751057743175429172241006516,
-0.0065136008537341835935751782245857190566213408663742,
-0.0065194204089507837716621608471744105332927968565350,
-0.0065229896301644093236200278475527046200258389534488,
-0.0065243147434844711101493981420344675346938325188668,
-0.0065234026697597116647835832553381463378208372112842,
-0.0065202610211247635437648564279952636659730883152448,
-0.0065148980973310139251876391822485156173881994046000,
-0.0065073228818632030269929887985300769493255117279122,
-0.0064975450378432519325984650757031567974217839709652,
-0.0064855749037228828907801581084588517065248075606706,
-0.0064714234887666620711349283646354271519108918189956,
-0.0064551024683271610864844978376034255243077846925030,
-0.0064366241789139993175879195166522620864125230016804,
-0.0064160016130585941723663220936775313635140314962246,
-0.0063932484139765108605939474097682212203675760921708,
-0.0063683788700293670449785451612908373679626674514068,
-0.0063414079089883108202809823514548926453557699541476,
-0.0063123510921011528533865852032236772154717900868154,
-0.0062812246079652951690632570670480581713043629367832,
-0.0062480452662086599688001730763456974955889980908304,
-0.0062128304909808820041509778354576785826700659869716};
	const mp_real k[]={0.096193022307830910744865275173058185724969796365840,
0.22080312441145242598386416451252109700887461915132,
0.34614911651644048867816794850643786742263180925102,
0.47166137756057126454972179647701835688091698798780,
0.59723670833951143791050375989554728831663400462532,
0.72284255871643690172591531902472706240995946988004,
0.84846546519517035836313573402105227344727235903904,
0.97409886122997210948231779052715628738277490700596,
1.0997391652816101918350915293842965818611827544220,
1.2253842587372790047019831570741693094909428651671,
1.3510328085429427473695418538685891209573082297027,
1.4766839341465617591907797225309181054594281034230,
1.6023370305851295717919722949597789164139184780796,
1.7279916685270692143009629091497368868350439168353,
1.8536475348664725607474315551645139669860059402104,
1.9793043958959126869441104613271308907880895101216,
2.1049620736445998411700514152156629320048133173075,
2.2306204302007991724673397109384732252550135902082,
2.3562793570432376853133762653846274234263235542302,
2.4819387676090867953141000105860380929526153506500,
2.6075985920082744176254413497004926490262627081911,
2.7332587731942719308396921535993751214739048489091,
2.8589192641437493130025229542451880221348757530316,
2.9845800257480735153528216187734221681208341513749,
3.1102410252155222015095748756493477306033424886455,
3.2359022348455051145508857396349221142749890387865,
3.3615636310775276063151855339199182051085035766962,
3.4872251937456461460504722027619443539468851978582,
3.6128869054884192022287106711051163565794448629031,
3.7385487512777909668476236617592364302929576025167,
3.8642107180398507511248646929571158751770224156568,
3.9898727943472238588111991600053831120415971815100,
4.1155349701677917838568013709025830515709880018784,
4.2411972366580646204070766876752979032654389288396,
4.3668595859922153420826390805099862647318910758084,
4.4925220112197963850997802448871332554376248801644,
4.6181845061466775851247102677611838638938632046240,
4.7438470652349012686251753790097563864102727546072,
4.8695096835180385186093745393336357804118961835392,
4.9951723565293178418103651346761657451418541476452,
5.1208350802403329631905423756919177560032574971172,
5.2464978510085565843158445996610775117814801763684,
5.3721606655322186439741158733640152998075894894236,
5.4978235208113711115113064593332866298005879005296,
5.6234864141141718619395561336024227862728775950532,
5.7491493429475893013580287201908194223323990764964,
5.8748123050318659837653081943301327241315509059320,
6.0004752982781902996352189189629182745147489680992,
6.1261383207691156995036650244061431517808580955108,
6.2518013707413409527819798109254418679355885728636,
6.3774644465705258529396415164000720645116786468576,
6.5031275467578671007623923585176733286605778861116,
6.6287906699182008346676170695281450554560430954192,
6.7544538147694330274982534980932718863851998701772,
6.8801169801231280086160172602930290414085338853020,
7.0057801648761097226149119143641378282397860070272,
7.1314433680029508272593589676500302994547615687072,
7.2571065885492420317684879834820294016901482060588,
7.3827698256255487256446164899166808353169757729048,
7.5084330784019743912030106288320336761714862263368,
7.6340963461032608928712319579838933092219356764516,
7.7597596280043647915976158690994500120940732114940,
7.8854229234264565894484801072289983181083323988932,
8.0110862317332964713344834988507643027561810074196,
8.1367495523279458468499528086127733243702875499080,
8.2624128846497789462182897971360698554084172957148,
8.3880762281717630078992978666320509224224722965500,
8.5137395823979793100298226270987181464542051364064,
8.6394029468613605262798255153256997758982056206640,
8.7650663211216226986075652391293888805582159345076,
8.8907297047633725736513809751788900555508146557936,
9.0163930973943731959405272798508000249717273525188,
9.1420564986439525320479136276091516757957192814888,
9.2677199081615415512389518477000415009391809104464,
9.3933833256153296407923209306625432065030165578628,
9.5190467506910265142223025329314303783397981749356,
9.6447101830907209005998788435964405691233772920616,
9.7703736225318273023522630471146385306259617343188,
9.8960370687461129939235176275787908488993893067056,
10.021700521478798218817854611320050052862740169862,
10.147363980487723240187712340175433294166276818370,
10.273027445542576520972709331454606593743001451308,
10.398690916424178862922591102767996220025292113173,
10.524354392923818827708160103053485081450933322889,
10.650017874842635204753638169360283384327589594184,
10.775681361991042685526935845534140014171814823202,
10.901344854188197258140573958773957819456394527063,
11.027008351261498153915911798757855749758816570673,
11.152671853046123463138470954579195500946570635649,
11.278335359384596794172668387071606459947463163362,
11.403998870126382581565223509721861401966715387432,
11.529662385127507857524695800828994099686250398972,
11.655325904250208489667562002768857879274930565696,
11.780989427362598058329327084436615175428742412664,
11.906652954338357700961967626994095195822563633172,
12.032316485056445390866758754052753211282664007975,
12.157980019400823244242311190159268371429047922276,
12.283643557260201564590951880013041995600957111456,
12.409307098527798438103843721744560896520004586500,
12.534970643101113788787609804076180508917658051856,
12.660634190881716888739100906970351591498252704621,
12.786297741775046397955493888710564750940359055886,
12.911961295690222080129901793947992908320214258431,
13.037624852539867406694819558230225321591255491652,
13.163288412239942321533720796321347852845073672666,
13.288951974709585493819596731933571218274303281966,
13.414615539870965436837736013808930536087920685619,
13.540279107649139916839179424924876206860365896298,
13.665942677971923118338090350276712897771885249007,
13.791606250769760071159014443455388415941824336764,
13.917269825975607880272327100004786196546733127689,
14.042933403524823332310878539630392138119460115464,
14.168596983355056482893149142458367182628045848723,
14.294260565406149856718650149119658982799274691162,
14.419924149620042918058277735827783774513081867251,
14.545587735940681492924476168626170160300209429138,
14.671251324313931846044241767609869205660918331373,
14.796914914687499135927097422486223001522306018951,
14.922578507010849989960733840459647250164287014937,
15.048242101235138958706685441972430300252394006120,
15.173905697313138624523195996647100005716211272050,
15.299569295199173154417861030392003716877365257186,
15.425232894849055100724783768235049544514686761113,
15.550896496220025265897382572532575054036702531844,
15.676560099270695459488456167409059178978179658811,
15.802223703960993986326482556960942233285406355125,
15.927887310252113715057912006991535765271258786424,
16.053550918106462585670200421017973175733117126654,
16.179214527487616423395110388503145908216336919167,
16.304878138360273934567263705847534198118714590051,
16.430541750690213767625661478164752864400900651913,
16.556205364444253529538605675970169735416354184776,
16.681868979590210654544330052581463285933102784354,
16.807532596096865028266628055435953468703265817700,
16.933196213933923276019865702849932793887657723500,
17.058859833071984629491328139269311241185544631153,
17.184523453482508291008769130716727426858294754432,
17.310187075137782219292984716340583909569738623843,
17.435850698010893264982852050782383564049800075788,
17.561514322075698588325338219836532413924614892334,
17.687177947306798295265565763764904701376988012938,
17.812841573679509231770638345112420008144508588920,
17.938505201169839879592671210411145267655005900168,
18.064168829754466299837118890795186767894638103727,
18.189832459410709073666629756531339907989132559753,
18.315496090116511192251765830329297864502355115073,
18.441159721850416850690478363641805322233472679426,
18.566823354591551103069768294966885846163436261283,
18.692486988319600338146175461624608755229890189091,
18.818150623014793537286535763037435608247452592102,
18.943814258657884278346003124779789313476459454402,
19.069477895230133451075169299250200844237669380529,
19.195141532713292651450132214172901122411912811049,
19.320805171089588224015900302601262703764818324880,
19.446468810341705922931392022566291480895114151452,
19.572132450452776163909831234857151987282447188711,
19.697796091406359840667440366522953300252788668345,
19.823459733186434680831473750838238297748184354272,
19.949123375777382117520912943902192493517876484333,
20.074787019163974654004312771201258169645692469923,
20.200450663331363699963765019136051175136149591861,
20.326114308265067858955858894271453394480252665955,
20.451777953950961647663707817352297708717054523642,
20.577441600375264628482166762508694833620023027060,
20.703105247524530937874629737355540605472521168421,
20.828768895385639193787397582900182831047100016150,
20.954432543945782766209460404319029931933367371905,
21.080096193192460395724373613999726479456206620480,
21.205759843113467145619271011790864947324187991431,
21.331423493696885673796336590397284045204605295882,
21.457087144931077811376479302584165756700914809474,
21.582750796804676435495609477779846278505732164018,
21.708414449306577624372756618165598546895326679043,
21.834078102425933083278126880987190404674531742261,
21.959741756152142830549790357789489908642615030336,
22.085405410474848133301621605195857050085181825974,
22.211069065383924682933899987952952822190780036320,
22.336732720869476001003024152023796954476563599686,
22.462396376921827066429435174451495705880063954561,
22.588060033531518155424322203015725220965673113643,
22.713723690689298885897173406821218195234767044862,
22.839387348386122458468833660896832769267290051743,
22.965051006613140086559472235492880409876024377906,
23.090714665361695608348720628896712549375045529550,
23.216378324623320273717126488071922049958435713192,
23.342041984389727699574844100775450050355291508232,
23.467705644652808987265954369733896399243035219797,
23.593369305404627996005739241626328643000936717510,
23.719032966637416766564344616453658704202335573152,
23.844696628343571089654227576914309029814147785190,
23.970360290515646213711235154328969226329536996398,
24.096023953146352686980703058947886469305198894065,
24.221687616228552329031159863623528053266401989670,
24.347351279755254327019608921949716732880607874054,
24.473014943719611452224440551777560640883148775532,
24.598678608114916392545276143355262215873808832317,
24.724342272934598196843912699522376148443501873384,
24.850005938172218827167444758101235253067369669152,
24.975669603821469815053991099299901117305612502788,
25.101333269876169018273624469037562973674458921332,
25.226996936330257474502451386151387739191626564241,
25.352660603177796348566654082131052646051888483859,
25.478324270412963970026007536239626204734089220644,
25.603987938030052957993223942108796729320830550363,
25.729651606023467430206741047051490822612776588707,
25.855315274387720293490530610624748936306109229567,
25.980978943117430612845415286193253098549121826622,
26.106642612207321056522489506466178787750367207787,
26.232306281652215414530772616412854380090869351060,
26.357969951447036188128398639390625679707369523768,
26.483633621586802247939673396374908905836619682146,
26.609297292066626558429402197426666439013101604025,
26.734960962881713966551195830113932372083303531082,
26.860624634027359052468175359082804559663385613421,
26.986288305498944040322784572717352416965965816639,
27.111951977291936767107441488205269574725781353265,
27.237615649401888707759667821779428100452599207017,
27.363279321824433054674270801843125372328070682919,
27.488942994555282849891251010994832975176221424564,
27.614606667590229168281502149823668362687332231088,
27.740270340925139350113176350898389847731151344970,
27.865934014555955281439928482749872814618900945540,
27.991597688478691720808235562658774757770604924213,
28.117261362689434670834718288854897119452150921233,
28.242925037184339793255970996607159002995994977403,
28.368588711959630866102929346134388842489702035123,
28.494252387011598281699362447353818651520139975885,
28.619916062336597584229754221420168479960760186261,
28.745579737931048045665719745923302892498102226859,
28.871243413791431278882264350310376783593230661425,
28.996907089914289886835710812603541136961136952551,
29.122570766296226146714064114114467962034684034663,
29.248234442933900728008021438512564915029692769555,
29.373898119824031443486831876919039713527922088519,
29.499561796963392032097827024660122768745223456630,
29.625225474348810972841738865908803454775168604804,
29.750889151977170328707950857881347184935582881099,
29.876552829845404619784645182229526694877490126233,
30.002216507950499724688464510769250960540441279226,
30.127880186289491809486848788750370685745568755490,
30.253543864859466283313682710724809128125584888563,
30.379207543657556779905341887007255863069264404627,
30.504871222680944164309697308116543931299085926716,
30.630534901926855564045168844340582680713545423452,
30.756198581392563424010547593346314490319080148644,
30.881862261075384584469070616275818670112751191827,
31.007525940972679381452165051445436238310638420965,
31.133189621081850768949415276964172935252654578074,
31.258853301400343462271678739743485387754834519030,
31.384516981925643101993913898184065023053305815944,
31.510180662655275437903216715777616649993214003646,
31.635844343586805532395818290239600653771843499781,
31.761508024717836982784402286417134341084940466044,
31.887171706046011161994082479078412684148460339045,
32.012835387569006477141762415609499780952389477258,
32.138499069284537645509404436185510700494064086444,
32.264162751190354987436986492724746319117480896896,
32.389826433284243735675643882257548298977414843665,
32.515490115564023360755699735867946915656582748263,
32.641153798027546911938002590170749825760644455203,
32.766817480672700373330230492007461557446581643204,
32.892481163497402034762606932998939283319694125380,
33.018144846499601877029821806009406523454357324586,
33.143808529677280971117877125039940130486048165944,
33.269472213028450891046098368274534237041918283218,
33.395135896551153139965683246904929716367557080720,
33.520799580243458589166915096813658041153626060626,
33.646463264103466929657561962140558567223142046483,
33.772126948129306135985028241060051754805270273722,
33.897790632319131941984536398077590913941118379766,
34.023454316671127328145004092907881972004132291846,
34.149118001183502020293359011538325423746557391890,
34.274781685854491999306811109387564212300050763199,
34.400445370682359021571090830483621417296958760186,
34.526109055665390149910872653220219282064304483070,
34.651772740801897294726546116445440797593906253626,
34.777436426090216765079180982532986615162479378049,
34.903100111528708829472968696335280162760763885179,
35.028763797115757286091617733193954609406315724040,
35.154427482849769042252144376924683699512477583080,
35.280091168729173702846241175113870928955852726130,
35.405754854752423167545930707474411938997827177930,
35.531418540917991236556529988474222345140610426471,
35.657082227224373224706068127835372576213433100756,
35.782745913670085583666223828898190595552919031854,
35.908409600253665532105586681456652936604797284828,
36.034073286973670693581603509146461314236881208750,
36.159736973828678741982954519721770333312788494638,
36.285400660817287054339319699112440426988848500944,
36.411064347938112370820549578011272977620878813010,
36.536728035189790461752151754810333922399938750923,
36.662391722570975801478750742818958319749886245086,
36.788055410080341248911778983217662443260614232112,
36.913719097716577734602116195258901064782055000544,
37.039382785478393954182717403574034897128131546412,
37.165046473364516068030461592993655605853156959567,
37.290710161373687407000517427285993034264970951206,
37.416373849504668184090464099136546381757193669620,
37.542037537756235211895228266801511066349830505790,
37.667701226127181625717606139530535364596752673483};

	
	//gettimeofday(&tpstart,NULL); 
	start = clock();

	char line[101];
	int nT;
	//ifstream data("Data.dat");
	ifstream data(argv[3]);
	ofstream result("result.dat");
	ofstream figure("figures.dat");
	ofstream bowls("bowls.dat");		//For deflection bowls
	
	cout<<"Read the file..."<<endl;
	data.getline(line,100);
	data>>nr_sys;
	data.getline(line,100);
	
    for(sys_nr=0; sys_nr<nr_sys; sys_nr++)
    {	
	result<<sys_nr<<endl;
	figure<<sys_nr<<endl;
	
	data.getline(line,100);
	data.getline(line,100);
	data.getline(line,100);
	data>>P0>>T0>>nT>>a>>layers>>nr_points;

	data.getline(line,100);
	data.getline(line,100);
	data.getline(line,100);

	points=new double[nr_points];
	E=new double[layers];
	nu=new double[layers];
	rho=new double[layers];
	h=new double[layers];
	friction=new int[layers];
	layerindex=new int[layers];

	for (i=0;i<nr_points;i++)
	{
		data>>points[i];
	}
	data.getline(line,100);
	data.getline(line,100);
	data.getline(line,100);
	data.getline(line,100);
	for (i=0;i<layers;i++)
	{
		data>>E[i]>>nu[i]>>rho[i]>>h[i]>>friction[i];
		data.getline(line,100);
		bowls<<E[i]<<"\t";	//For deflection bowls
	}
	//data.close();
	cout<<(sys_nr+1)<<"\t"<<"start!"<<endl;

	/* layer matrix index*/
	layerindex[0]=0;
	for (i=1;i<layers;i++)
	{
		if (friction[i-1]==1) layerindex[i]=layerindex[i-1]+2;
		else layerindex[i]=layerindex[i-1]+3;
	}	
	nr_layerindex=layerindex[layers-1]+2;
	/***********************/

	cout<<"Assign value for the variables...";
	R=25.0;
	Tt=nT*T0;
	t0=2.0*T0;
	t=new double[N];
	
	for (i=0;i<N;i++)
	{
		t[i]=(i+1)*Tt/N;
		if ((t[i]>t0) && (t[i]<(t0+T0)))
			Ft[i]=sin((t[i]-t0)/T0*Pi);
		else
			Ft[i]=mp_complex(0.0,0.0);
	}
	
	for (i=0;i<N;i++)
	{
		Fn[i]=mp_complex(0.0,0.0);
		for (j=0;j<N;j++)
			Fn[i]+=Ft[j]*exp(-2*Pi*(i+1)*(j+1)/N*(I));
		Fn[i]/=N;
	}
	
	for (i=0;i<N;i++) omega[i]=2*Pi*(i+1)/Tt;
	
	mu=new mp_real[layers];
	lambda=new mp_real[layers];
	cp=new mp_real[layers];
	cs=new mp_real[layers];
	G=new mp_real[layers];
	for (layers_nr=0;layers_nr<layers;layers_nr++)
	{
		mu[layers_nr]=E[layers_nr]/2/(1+nu[layers_nr]);
 		lambda[layers_nr]=nu[layers_nr]*E[layers_nr]/(1+nu[layers_nr])/(1-2*nu[layers_nr]);
 		cp[layers_nr]=sqrt((lambda[layers_nr]+2*mu[layers_nr])/rho[layers_nr]);
 		cs[layers_nr]=sqrt(mu[layers_nr]/rho[layers_nr]);
 		G[layers_nr]=E[layers_nr]/(2*(1+nu[layers_nr]));
	}
	
	P=new mp_complex[nr_layerindex];//*************************
	Ksys = new mp_complex *[nr_layerindex];//************************
    for (i=0;i<(nr_layerindex);i++)
    {
      	Ksys[i]=new mp_complex [nr_layerindex];
    }
	w=new mp_complex **[N];
	for (i=0;i<N;i++)
	{
		w[i]=new mp_complex *[M];
		for (j=0;j<M;j++)
		{
			w[i][j]=new mp_complex [2*layers];
		}
	}
	for (i=0;i<nr_layerindex;i++)//******************************
	{
		for (j=0;j<nr_layerindex;j++)
		{
			Ksys[i][j]=mp_complex(0.0,0.0);
		}
	}
	cout<<"done!"<<endl;

	//Assemble
	cout<<"Assembling...";
	for (n=0;n<N;n++)
	{
		for (m=0;m<M;m++)
		{
			for (i=0;i<nr_layerindex;i++)//*******************************
			{
				for (j=i;j<nr_layerindex;j++)
					Ksys[i][j]=mp_complex(0.0,0.0);
			}
			if (layers>1)
			{
				for (layers_nr=0;layers_nr<layers-1;layers_nr++)
				{
					
					s=mp_complex(0.0,0.0)-sqrt(mp_complex(1.0,0.0)-pow(omega[n],2)/pow(cp[layers_nr]*k[m],2))*I;
 					ts=mp_complex(0.0,0.0)-sqrt(mp_complex(1.0,0.0)-pow(omega[n],2)/pow(cs[layers_nr]*k[m],2))*I;
					sin_s=sin(k[m]*s*h[layers_nr]);
					cos_s=cos(k[m]*s*h[layers_nr]);
					sin_ts=sin(k[m]*ts*h[layers_nr]);
					cos_ts=cos(k[m]*ts*h[layers_nr]);
					DD=mp_complex(2.0,0.0)*(mp_complex(1.0,0.0)-cos_s*cos_ts)+(s*ts+mp_complex(1.0,0.0)/(s*ts))*sin_s*sin_ts;

					K[0][0]=((mp_complex(1.0,0.0)+pow(ts,2))*k[m]*G[layers_nr])/DD*(mp_complex(1.0,0.0)/ts*cos_s*sin_ts+s*sin_s*cos_ts);
					K[0][1]=((mp_complex(1,0)+pow(ts,2))*k[m]*G[layers_nr])/DD*((mp_complex(3,0)-pow(ts,2))/(mp_complex(1,0)+pow(ts,2))*(mp_complex(1,0)-cos_s*cos_ts)+(mp_complex(1,0)+mp_complex(2,0)*pow(s,2)*pow(ts,2)-pow(ts,2))/(s*ts*(mp_complex(1,0)+pow(ts,2)))*sin_s*sin_ts);
					K[0][2]=((mp_complex(1.0,0.0)+pow(ts,2))*k[m]*G[layers_nr])/DD*(-s*sin_s-mp_complex(1.0,0.0)/ts*sin_ts);
					K[0][3]=((mp_complex(1.0,0.0)+pow(ts,2))*k[m]*G[layers_nr])/DD*(cos_s-cos_ts);
					K[1][1]=((mp_complex(1.0,0.0)+pow(ts,2))*k[m]*G[layers_nr])/DD*(mp_complex(1.0,0.0)/s*sin_s*cos_ts+ts*cos_s*sin_ts);
 					K[1][3]=((mp_complex(1.0,0.0)+pow(ts,2))*k[m]*G[layers_nr])/DD*(mp_complex(-1.0,0.0)/s*sin_s-ts*sin_ts);
 					K[1][2]=-K[0][3];
 					K[2][2]=K[0][0];
 					K[2][3]=-K[0][1];
 					K[3][3]=K[1][1];
					/*for(i=0;i<4;i++)
					{
						for(j=i;j<4;j++)
						{
							Ksys[2*(layers_nr)+i][2*(layers_nr)+j] += K[i][j];
						}
					}*/
					/*new*/
					if (friction[layers_nr]==1)
					{
						for(i=0;i<4;i++)
						{
							for(j=i;j<4;j++)
							{
								Ksys[layerindex[layers_nr]+i][layerindex[layers_nr]+j] += K[i][j];
							}
						}
					}
					else
					{
						for(i=0;i<3;i++)
						{
							for(j=i;j<3;j++)
							{
								Ksys[layerindex[layers_nr]+i][layerindex[layers_nr]+j] += K[i][j];
							}
							Ksys[layerindex[layers_nr]+i][layerindex[layers_nr]+4]+=K[i][3];
						}
						Ksys[layerindex[layers_nr]+4][layerindex[layers_nr]+4]+=K[3][3];
					}
				}
			}
			s0=mp_complex(0.0,0.0)-I*sqrt(mp_complex(1.0,0.0)-pow(omega[n],2)/pow(cp[layers-1]*k[m],2));
			ts0=mp_complex(0.0,0.0)-I*sqrt(mp_complex(1.0,0.0)-pow(omega[n],2)/pow(cs[layers-1]*k[m],2));
			K0[0][0]=k[m]*G[layers-1]*I*s0*(mp_complex(1.0,0.0)+pow(ts0,2))/(mp_complex(1.0,0.0)+s0*ts0);
			K0[0][1]=k[m]*G[layers-1]*2-(mp_complex(1.0,0.0)+pow(ts0,2))/(mp_complex(1.0,0.0)+s0*ts0);
			K0[1][1]=k[m]*G[layers-1]*I*ts0*(mp_complex(1.0,0.0)+pow(ts0,2))/(mp_complex(1.0,0.0)+s0*ts0);
			for(i=0;i<2;i++)
			{
				for(j=i;j<2;j++)
					Ksys[nr_layerindex-2+i][nr_layerindex-2+j] += K0[i][j];
			}
			//symmetric
			for (i=0;i<nr_layerindex;i++)//**************************
			{
				for (j=i-1;j>=0;j--)
				{
					Ksys[i][j]=Ksys[j][i];
				}
			}
			for (i=0;i<nr_layerindex;i++) P[i]=mp_complex(0.0,0.0);
			P[1]=mp_complex(1.0,0.0);
			GSolve(Ksys,nr_layerindex,P);
			/*save all displacement data, for stress calculating later if necessary,
			which is not included in the FWD simulation*/
			for (i=0;i<layers*2;i++) w[n][m][i]=P[i];
		}
	}
	cout<<"done!"<<endl;

	//Inverse Transform
	cout<<"Start Inverse Transforming...";	
	Result=new mp_real *[nr_points];
	for (i=0;i<nr_points;i++) Result[i]=new mp_real[2];
	Figure=new mp_real *[N];
	for (i=0;i<N;i++) Figure[i]=new mp_real[nr_points+2];
		
	for(n=0;n<N;n++)
	{
		Figure[n][0]=t[n];
 		Figure[n][1]=(P0*Ft[n]).real;
	}

	for(points_nr=0;points_nr<nr_points;points_nr++)
	{
		r=points[points_nr]/100;
		for(i=0;i<M;i++) T[i]=bessj0(k[i]*r);
		for(n=0;n<N;n++)
		{
			W[n]=mp_complex(0.0,0.0);			
			for(m=0;m<M;m++)
			{
				/*surface deflections
				w[n][m][0]--horizontal displacement at the top of FIRST layer
				w[n][m][1]--vertical displacement at the top of FIRST layer
				w[n][m][2]--horizontal displacement at the bottom of FIRST layer
				w[n][m][3]--vertital displacement at the bottom of FIRST layer
				w[n][m][4]--horizontal displacement at the bottom of SECOND layer, if full bonding between layers
					  --horizontal displacement at the top of SECOND layer, if no bonding between layers
				vertical displacements at bottem of first layer and at top of the second layer are identical
				*/
				W[n]+=Fm[m]*T[m]*w[n][m][1];
			}
			W[n]*=Fn[n];
		}

		for(i=0;i<N;i++)
		{
			Wt[i]=mp_complex(0.0,0.0);
			for(j=0;j<N;j++)
			{
				Wt[i] += W[j]*exp(2*Pi*(i+1)*(j+1)/N*(I));
			}
			Wt[i] *= mp_complex(-P0,0.0);
		}
		for(i=0;i<N;i++) Figure[i][points_nr+2]=Wt[i].real;
		Min=0;
		time_Min=0;
		for(i=0;i<N;i++)
		{
			if(Wt[i].real<Min) 
			{
				Min=Wt[i].real;
				time_Min=t[i];
			}
		}
		Result[points_nr][0]=time_Min;
  		Result[points_nr][1]=Min;
		result<<Result[points_nr][0]<<"\t"<<Result[points_nr][1]<<endl;
		bowls<<Result[points_nr][1]<<"\t";
	}
	//result.close();
	result<<endl<<endl;
	bowls<<endl;
	for(n=0;n<N;n++)
	{
		for(points_nr=0;points_nr<nr_points+2;points_nr++) figure<<Figure[n][points_nr]<<"\t";
		figure<<endl;
	}
	//figure.close();
	figure<<endl<<endl;
	//cout<<"done!"<<endl;
	cout<<"done!"<<endl;

	/*gettimeofday(&tpend,NULL); 
	timeuse=1000000*(tpend.tv_sec-tpstart.tv_sec)+
	tpend.tv_usec-tpstart.tv_usec;
	timeuse/=1000000; 
	printf("Used Time:%f\n",timeuse); */
	
	/*stop = clock();//结束时间
	cout<<"Used Time:"<<(float)((stop-start)/1000)<<"Sencods"<<endl;
	cout << "Complete!" <<endl;
	cout<<"Press any key...";*/
	//getchar();
	
	//内存释放
	/*
	delete points;
	points=NULL;
	
	delete E;
	E=NULL;
	
	delete nu;
	nu=NULL;
	
	delete rho;
	rho=NULL;
	
	delete h;
	h=NULL;
	
	delete friction;
	friction=NULL;
	
	delete layerindex;
	layerindex=NULL;
	
	delete t;
	t=NULL;
	
	delete mu;
	mu=NULL;
	
	delete lambda;
	lambda=NULL;
	
	delete cp;
	cp=NULL;
	
	delete cs;
	cs=NULL;
	
	delete G;
	G=NULL;
	
	delete P;
	P=NULL;
	
    	for (i=0;i<(nr_layerindex);i++)
    	{
      		delete Ksys[i];
		Ksys[i]=NULL;
    	}
    	delete Ksys;
	Ksys=NULL;
	
	for (i=0;i<N;i++)
	{
		for (j=0;j<M;j++)
		{
			delete w[i][j];
			w[i][j]=NULL;
		}
		delete w[i];
		w[i]=NULL;
	}
	delete w;
	w=NULL;
	
	for (i=0;i<nr_points;i++) 
	{
		delete Result[i];
		Result[i]=NULL;
	}
	delete Result;
	Result=NULL;
	
	for (i=0;i<N;i++) 
	{
		delete Figure[i];
		Figure[i]=NULL;
	}
	delete Figure;
	Figure=NULL;
	*/
	delete [] points;	
	delete [] E;	
	delete [] nu;	
	delete [] rho;	
	delete [] h;	
	delete [] friction;	
	delete [] layerindex;	
	delete [] t;	
	delete [] mu;	
	delete [] lambda;	
	delete [] cp;	
	delete [] cs;	
	delete [] G;	
	delete [] P;
	
    	for (i=0;i<(nr_layerindex);i++)
    	{
      		delete [] Ksys[i];
    	}
    	delete [] Ksys;
	Ksys=NULL;
	
	for (i=0;i<N;i++)
	{
		for (j=0;j<M;j++)
		{
			delete [] w[i][j];
		}
		delete [] w[i];
	}
	delete [] w;
	w=NULL;
	
	for (i=0;i<nr_points;i++) 
	{
		delete [] Result[i];
	}
	delete [] Result;
	Result=NULL;
	
	for (i=0;i<N;i++) 
	{
		delete [] Figure[i];
	}
	delete [] Figure;
	Figure=NULL;
	
    }
	cout<<"Finished!"<<endl;
	result.close();
	figure.close();
	bowls.close();
	stop = clock();//结束时间
	cout<<"Used Time:"<<(float)((stop-start)/1000)<<"Sencods"<<endl;
	cout << "Complete!" <<endl;
	cout<<"Press any key...";
	mp::mp_finalize();
	return 0;
}