c CODE_TO_DONES/mcnewd_tilt.for c**************************************************************** program mcnewd_tilt c**************************************************************** c driver for calculating of photometric properties c uses mcnewd_sub.for for the calculation of scattering c from a single position-file / sun direction) c several observing directions corresponding to c a fixed phase-angle and several sun/observer latitudes c output of results via tiltfile (unit21) c use mc_tiltread_f.pro for reading the output c**************************************************************** c HS 7/12/99 - 22/01/02 implicit double precision (a-h,o-z) include 'mcnew_main.inc' parameter(nnt=1000) character*60 posfile,parfile,outfile,epsfile,wwwfile double precision testfii(nnt),testlat(nnt) double precision testw(nnt),testd(nnt),testn(nnt) double precision testm(nnt,10) double precision lonsun,mu0 double precision iscat_mean,npot_mean,nbox_mean character*8 date character*10 time character*5 zone integer*4 values(8) character*12 smet character*60 tiltfile double precision bsuntable(100) double precision testfii_ori(nnt),testlat_ori(nnt) ar=0.017453293d0 ra=57.29577951d0 c posfile = posdat-file containing particle positions c parfile = parameter-file (junk -> defaults) c outfile = output of weights + directions c epsfile = output of weights + fiilat c wwwfile = output of weights c sredu = reduction factor for particle radius (sredu < 1) c zredu = reduction factor for vertical coordinates (zredu > 1) c idum = seed c tiltfile = output of multiple weights c******************************************************** c READ INPUT VALUES c******************************************************** c particle positions from POSFILE c BSUN,LONSUN = elevation and longitude of sun c TESTLAT,TESTFII = table of observer positions c OUTPUT: c testw,testn,testd,testm c to convert testw to I/F: c multiply by sin(bsun)/nphot/sin(bobs) read(*,*) posfile read(*,*) bsun,lonsun read(*,*) ntest c testlat and testfii intepreted later do 10 i=1,ntest read(*,*) testlat_ori(i),testfii_ori(i) 10 continue read(*,*) nphot,mx,my read(*,*) parfile read(*,*) outfile read(*,*) epsfile read(*,*) wwwfile read(*,*) method smet='unknown' if(method.eq.1 ) smet='LAMBERT-element ' if(method.eq.11) smet='LAMBERT-sphere-center ' if(method.eq.12) smet='LAMBERT-sphere-point ' if(method.eq.2 ) smet='ISOTROPIC-sphere' if(method.eq.3 ) smet='HG-sphere-center' if(method.eq.4 ) smet='HG-sphere-point' if(method.eq.13) smet='HG2-sphere-center' if(method.eq.14) smet='HG2-sphere-point' if(method.eq.5 ) smet='EULER-sphere-center' if(method.eq.6 ) smet='EULER-sphere-point' if(method.eq.7 ) smet='CALLISTO-sphere-center' if(method.eq.8 ) smet='CALLISTO-sphere-point' if(method.eq.9 ) smet='POWER-sphere-center' if(method.eq.10) smet='POWER-sphere-point' read(*,*) ppar1,ppar2,ppar3 c icheck=-1 -> no check output read(*,*) icheck read(*,*) idum read(*,*) sredu,zredu read(*,*) sizemin read(*,*) s20,f20 c----------------------------------------- c Saturn-shine: c iglobe=0 -> SUN c iglobe=1 -> SATURN c iglobe=2 -> SATURN-SOUTH c iglobe=3 -> SATURN-NORTH c abox = calculation distance c rsaturn = Saturn radius c nlat_sat = number of latitude bins c nfii_sat = number of longitude bins c model_sat = saturn model c = 1 -> Lambert c = 2 -> RED (DONES ET ALL) c = 3 -> BLUE (DONES ET ALL) rsaturn=60.4 read(*,*) iglobe read(*,*) abox read(*,*) nlat_sat,nfii_sat read(*,*) model_sat if(model_sat.eq.1) smodel='Lambert-planet' if(model_sat.eq.2) smodel='RED 0.5*zones+0.5*belts' if(model_sat.eq.3) smodel='BLUE 0.5*zones+0.5*belts' if(model_sat.eq.4) smodel='RED+BLUE 0.5*zones+0.5*belts' if(model_sat.eq.5) smodel='ZONES 0.75*blue+0.25*red' if(model_sat.eq.6) smodel='BELTS 0.75*blue+0.25*red' c -- INTEPRET TILT INPUT VALUES c multiple values of sunlat and obslat c alpha_choice = 0 -> testfii intepreted as dfii c 1 -> testfii intepreted as alpha c lat_choice = 0 -> testlat intepreted as they are c lat_choice = 1 -> testlat added to bsun read(*,*) alpha_choice,lat_choice read(*,*) nlatsun do 3 i=1,nlatsun read(*,*) bsuntable(i) 3 continue read(*,*) tiltfile write(*,*) '==========================================' write(*,*) 'CODE_TO_DONES/mcnewd_tilt.for (22.01.2002)' call date_and_time (date, time, zone, values) write(*,999) values(1),values(2),values(3),values(5),values(6) 999 format(1x,i4,'-',i2,'-',i2,2x,i2,':',i2,':',i2) write(*,*) '==========================================' write(*,*) 'POSFILE= ',posfile write(*,*) 'BSUN, LONSUN=',bsun,lonsun write(*,*) 'NTEST ',ntest write(*,*) 'NPHOT,MX,MY= ',nphot,mx,my write(*,*) 'METHOD = ',method,' ',smet write(*,*) 'PPAR 1,2,3 = ',ppar1,ppar2,ppar3 write(*,*) 'SEED = ',idum write(*,*) 'parfile,outfile,epsfile,wwwfile=' write(*,*) ' -> ',parfile write(*,*) ' -> ',outfile write(*,*) ' -> ',epsfile write(*,*) ' -> ',wwwfile write(*,*) 'sredu,zredu',sredu,zredu write(*,*) 'sizemin ',sizemin write(*,*) 's20,f20 ',s20,f20 if(iglobe.eq.0) write(*,*) 'ILLUMINATION BY SUN' if(iglobe.eq.1) write(*,*) 'ILLUMINATION BY SATURN' if(iglobe.eq.2) write(*,*) 'ILLUMINATION BY SATURN-SOUTH' if(iglobe.eq.3) write(*,*) 'ILLUMINATION BY SATURN-NORTH' if(iglobe.ge.1) then write(*,*) 'MODEL_SAT=',model_sat,smodel write(*,*) 'DISTANCE =',abox,abox/rsaturn write(*,*) 'NLAT,NFII=',nlat_sat,nfii_sat endif write(*,*) 'alpha_choice=',alpha_choice write(*,*) 'lat_choice =',lat_choice write(*,*) 'alpha_choice = 0 -> testfii intepreted as dfii' write(*,*) 'alpha_choice = 1 -> testfii intepreted as alpha' write(*,*) 'lat_choice = 0 -> testlat as it is' write(*,*) 'lat_choice = 1 -> testlat added to bsun' write(*,*) 'nlatsun =',nlatsun write(*,*) (bsuntable(i),i=1,nlatsun) write(*,*) tiltfile write(*,*) '==========================================' open(unit=21,file=tiltfile,status='unknown') versio21=1 write(21,*) versio21 write(21,*) posfile write(21,*) sunfii,sunlat,nphot write(21,*) method write(21,*) ppar1,ppar2,ppar3 write(21,*) seed write(21,*) sredu,zredu write(21,*) sizemin write(21,*) s20,f20 write(21,*) ntest write(21,*) (testfii_ori(i),i=1,ntest) write(21,*) (testlat_ori(i),i=1,ntest) write(21,*) model_sat write(21,*) iglobe write(21,*) abox,abox/rsaturn write(21,*) nlat_sat,nfii_sat write(21,*) alpha_choice,lat_choice write(21,*) nlatsun write(21,*) (bsuntable(i),i=1,nlatsun) write(21,*) tiltfile c multiple values of sunlat and obslat c alpha_choice = 0 -> testfii intepreted as dfii c 1 -> testfii intepreted as alpha c lat_choice = 0 -> testlat intepreted as they are c lat_choice = 1 -> testlat added to bsun c cos_alpha=cos(dlobs*ar)*cos(bsun*ar)*cos(bobs*ar)+ c sin(ar*bobs)*sin(bsun*ar) do 34 ilatsun=1,nlatsun sunlat=bsuntable(ilatsun) sunfii=lonsun write(21,*) sunlat,sunfii write(*,*) '-----------------------------' write(*,*) 'sun: lat,fii=',sunlat,sunfii do 35 itest=1,ntest if(lat_choice.eq.0) + testlat(itest)=testlat_ori(itest) if(lat_choice.eq.1) + testlat(itest)=testlat_ori(itest)+sunlat if(alpha_choice.eq.0) testfii(itest)=testfii_ori(itest) if(alpha_choice.eq.1) then apu1=cos(testfii_ori(itest)*ar)- + sin(ar*sunlat)*sin(ar*testlat(itest)) apu2=cos(ar*sunlat)*cos(ar*testlat(itest)) apu3=apu1/apu2 if(apu3.ge.1) apu3=.999999 if(apu3.le.-1) apu3=-.999999 apu3=acos(apu3)*ra if(testfii_ori(itest).eq.0) apu3=0. testfii(itest)=sunfii+apu3 endif apu1=sin(ar*sunlat)*sin(ar*testlat(itest)) apu2=cos(ar*sunlat)*cos(ar*testlat(itest))* + cos((sunfii-testfii(itest))*ar) apu=apu1+apu2 if(apu.gt.0.999999) apu=0.999999 if(apu.lt.-0.999999) apu=-0.999999 alpha=acos(apu)*ra write(21,*) testlat(itest),testfii(itest),alpha write(*,*) 'obs: lat,fii,alpha ', + testlat(itest),testfii(itest),alpha 35 continue c Saturn-shine ? nphot_tod=nphot flux_sat=0. if(iglobe.ge.1) then call saturn_illumination + (sunlat,sunfii,nphot_ori,nphot_tod,flux_sat) nphot=nphot_tod endif write(*,*) 'flux_sat,nphot',flux_sat,nphot c MC-calculation call tim(cpu1) call mcnew_sub(nphot,posfile,parfile,outfile, + sunfii,sunlat, + ntest,testfii,testlat,testw,testd,testn,testm, + idum,icheck,mx,my, + npart,npartori,npot_mean,nbox_mean,iscat_mean, + method,epsfile,wwwfile,ifor,sredu,tau,ndirect, + zredu,ppar1,ppar2,ppar3) call tim(cpu2) mu0=cos((90.d0-sunlat)*0.017453293) fdirect=1.*ndirect/nphot tau_path=-log(fdirect) tau_n=tau_path*mu0 cpu_comp=(cpu2-cpu1)/(1.*nphot/1000.) write(*,*) 'NORI,TAU = ',npartori,tau write(*,*) 'SREDU, ZREDU = ', sredu,zredu write(*,*) 'cpu, cpu_comp',cpu2-cpu1,cpu_comp write(*,*) 'npart,npot,nbox,iscat', + npart,npot_mean,nbox_mean,iscat_mean write(*,*) 'fdir,mu0.tau_p,tau_n', + fdirect,mu0,tau_path,tau_n write(*,*) 'sizemin ',sizemin write(*,*) 's20,f20 ',s20,f20 write(21,*) (testw(i),i=1,ntest) write(21,*) (testn(i),i=1,ntest) write(21,*) (testd(i),i=1,ntest) write(21,*) ((testm(i,j),i=1,ntest),j=1,10) write(21,*) iscat_mean,npot_mean,nbox_mean write(21,*) npartori,npart,mx,my write(21,*) tau write(21,*) fdirect,tau_path,tau_n write(21,*) cpu2-cpu1,cpu_comp 34 continue close(21) write(*,*) '==========================================' call date_and_time (date, time, zone, values) write(*,999) values(1),values(2),values(3),values(5),values(6) write(*,*) '==========================================' end